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Chapter  1 

Introduction 


The  project  entitled  ^"Large  Deformation  Analysis  of  Nonlinear  Homogeneous  and  Heterogeneous 
Media  Using  an  Adaptive  Arbitrary  Lagrangian-Eulerian  Finite  Element  Method^  at  the  University 
of  Alabama  in  March  1991,  where  the  PI  was  then  an  Assistant  Professor  of  Engineering  Mechan¬ 
ics.  It  was  subsequently  transferred  to  the  Ohio  State  University,  where  the  PI  relocated  beginning 
September  1991.  A  no-cost  extension  was  granted  to  continue  this  project  till  September  30,  1995. 
During  the  four  year  period  of  this  grant,  substantial  progress  has  been  made  in  advancing  the 
state  of  the  art  in  multiple  scale  modeling  of  advanced  heterogeneous  materials,  in  large  defor¬ 
mation  analysis  of  solids  with  applications  in  metal  forming,  and  also  in  solidification  modeling. 
Research  has  been  conducted  in  a  few  distinct  areas  that  are  delineated  below.  A  list  of  publica¬ 
tions,  acknowledging  this  grant  is  provided  in  chapter  2. 

I.  The  Voronoi  Cell  Finite  Element  Model  for  Heterogeneous  Microstructures 
( Details  in:  Chapter  3) 

A  new  finite  element  model  has  been  developed  for  analysis  of  heterogeneous  media,  in  which 
the  second  phase  is  randomly  dispersed  within  the  matrix,  A  mesh  generator  based  on  Dirichlet 
tessellation  discretizes  the  domain,  accounting  for  arbitrariness  in  location,  shape  and  size  of  the 
second  phase.  This  results  in  a  network  of  convex  “Voronoi  polygons”,  which  forms  the  elements  in 
a  finite  element  mesh.  An  assumed  stress  hybrid  formulation  has  been  implemented  for  accommo¬ 
dating  arbitrary  multi-sided  elements  in  the  finite  element  model.  Composite  element  formulations 
have  been  developed  to  incorporate  the  effect  of  second  phase  within  each  element.  Formulations 
have  been  developed  for  elasticity,  thermo-elasticity,  elasto-plasticity  and  heat  conduction. 

II.  Multiple  Scale  Modeling  of  Heterogeneous  Materials  Using  Asymptotic  Homoge¬ 
nization  and  the  Voronoi  Cell  FEM 

(Details  in:  Chapter  4) 

In  this  work,  a  multiple  scale  finite  element  model  (VCFEM-HOMO)  has  been  developed  for 
elastic-plastic  analysis  of  heterogeneous  (porous  and  composite)  materials  by  combining  asymp¬ 
totic  homogenization  theory  with  the  Voronoi  CeU  finite  element  model  (VCFEM).  VCFEM  for 
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microstructural  modeling  originates  from  Dirichlet  tessellation  of  representative  material  elements 
at  sampling  points  in  the  structure.  Structural  modeling  is  done  by  the  general  purpose  finite  ele¬ 
ment  code  ABAQUS,  and  interfacing  with  the  microscale  VCFExM  analysis  is  done  through  the  user 
subroutine  in  ABAQUS  for  material  constitutive  relation,  LEMAT.  Asymptotic  homogenization  in 
UMAX  generates  macroscopic  material  parameters  for  ABAQUS,  Following  the  macroscopic  analy¬ 
sis,  a  local  VCFEM  analysis  is  invoked  to  depict  the  true  evolution  of  microstructural  state  variables. 
Various  numerical  examples  are  executed  for  validating  the  effectiveness  of  VCFEM-HOMO,  and 
the  effect  of  size,  shape  and  distribution  of  heterogeneities  on  local  and  global  response  is  examined. 

III.  R-S  Adapted  Arbitrary  Lagrangian-Eulerian  Finite  Element  Method  For  Metal 
Forming  Problems  with  Strain  Localization 

(Details  in:  Chapter  5) 

In  this  work,  an  adaptive  arbitrary  Lagrangian-Eulerian  (ALE)  finite  element  method  is  devel¬ 
oped  for  solving  large  deformation  problems,  with  applications  in  metal  forming  simulation.  The 
ALE  mesh  movement  is  coupled  with  r-adaptation  of  automatic  node  relocation,  to  minimize  ele¬ 
ment  distortion  during  the  process  of  deformation.  Strain  localization  is  considered  in  this  study 
through  the  constitutive  relations  for  ductile  porous  materials,  proposed  by  Gurson  and  Tvergaard. 
Prediction  of  localized  deformation  is  achieved  through  a  multi-level  mesh  superimposition  method, 
termed  as  s-adaptation.  The  model  is  validated  by  comparison  with  established  results  and  codes, 
and  a  few  metal  forming  problems  are  simulated  to  test  its  effectiveness. 


IV,  Arbitrary  Lagrangian-Eulerian  Finite  Element  Method  in  Solidification  Model¬ 
ing 

( Details  in:  Chapter  6) 

A  heat  transfer  analysis  for  solidification  problems  has  been  conducted  to  evaluate  the  temperature- 
field  and  the  location  of  the  phase  change  interface.  The  arbitrary  Lagrangian-Eulerian  kinematic 
description  has  been  utilized  in  the  finite  element  formulation  for  imparting  flexibility  to  the  motion 
of  the  nodes  in  model.  By  detaching  the  nodal  points  from  the  underlying  material,  nodes  can  be 
monitored  to  follow  the  evolving  front  while  maintaining  shapes  of  the  elements.  Special  numeri¬ 
cal  techniques  to  smoothen  the  deforming  front  and  to  avoid  continuous  remeshing  are  introduced. 
Numerical  examples  have  been  solved  to  establish  the  validity  of  the  present  model  and  its  strength. 
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Chapter  2 

List  of  Publications 


2.1  Theses  Acknowledging  this  Grant 

1.  “Voronoi  Cell  Finite  Element  Method  for  Micropolar  Thermoelastic  Heterogeneous  Materi¬ 
als”,  Y.S.  Liu,  M.S.  Thesis,  The  Ohio  State  University,  1994. 

2.  “R-S  Adapted  Arbitrary  Lagrangian  Eulerian  Finite  Element  Method  For  Metal  Forming 
Analysis  With  Strain  Localization”,  S.  Raju,  M.S.  Thesis,  The  Ohio  State  University,  1995. 

2.2  Publications  Acknowledging  this  Grant 

I,  Voronoi  Cell  Finite  Element  Model  for  Heterogeneous  Microstructures 

L  S.  Ghosh  and  S.  Moorthy,  “A  Model  for  Analysis  of  Arbitrary  Composite  and  Porous  Mi¬ 
crostructures  with  Voronoi  Cell  Finite  Elements,”  (in  review). 

2.  S.  Moorthy  and  S.  Ghosh,  ‘‘  Mesoscopic  Analysis  of  Small  Deformation  in  Heterogeneous 
Materials  using  Voronoi  Cell  Finite  Element  Model,”  Computational  Mechanics'  95  Theory  and 
Applications,  S.  N.  Atluri,  G.  Yagawa,  and  T.  Cruse  eds.,  pp  1916-1921,  1995. 

3.  S.  Ghosh  and  Y.  Liu,  “Voronoi  CeU  Finite  Element  Model  Based  on  Micropolar  Theory  of 
Thermoelasticity  for  Heterogeneous  Materials,”  International  Journal  for  Numerical  Methods  in 
Engineering,  Vol.  38,  pp  1361-1398,  1995. 

4.  S.  Ghosh  and  S.  Moorthy,  “Elastic-Plastic  Analysis  of  Heterogeneous  Microstructures  Us¬ 
ing  the  Voronoi  CeU  Finite  Element  Method,”  Computer  Methods  in  Applied  Mechanics  and 
Engineering,  Vol.  121,  pp  373-409,  1995. 

5.  S.  Ghosh  and  R.  L.  MaUett,  “  Voronoi  CeU  Finite  Elements,”  Computers  and  Structures,  Vol. 
50,  No.l,  pp  33-46,  1994. 

6.  S.  Moorthy,  S.  Ghosh  and  Y.S.  Liu,  ,  “  Voronoi  CeU  Finite  Element  Model  for  Thermo- 
Elastoplastic  Deformation  in  Random  Heterogeneous  Media,”  Applied  Mechanics  Reviews, 
Vol.  47,  No.  1,  Part  2,  pp  207-221,  1994. 
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7.  S.  Moorthy,  S.  Ghosh  and  Y.S.  Liu,  “  Voronoi  Cell  Finite  Element  Model  for  Random 
Micropolar  Elastic- Plastic  Heterogeneous  Media,”  Proceedings  of  13th  Army  Symposium  on 
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8.  S.  Ghosh  and  S.  Moorthy,“  A  Voronoi  Cell  Finite  Element  Model  for  Random  Heteroge¬ 
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9.  S.  Ghosh  and  S.  N.  Mukhopadhyay,  “  A  Material  Based  Finite  Element  Analysis  of  Het¬ 
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and  Engineering,  Vol.  104,  pp  211-247,  1993. 

II.  Multiple  Scale  Modeling  of  Heterogeneous  Materials  Using  Asymptotic 
Homogenization  and  the  Voronoi  Cell  FEM 


10.  S.  Ghosh,  K.  Lee  and  S.  Moorthy,  “Two  Scale  Analysis  of  Heterogeneous  Elastic-Plastic 
Materials  with  Asymptotic  Homogenization  and  Voronoi  Cell  Finite  Element  Model,”  (in 
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Chapter  3 

The  Voronoi  Cell  Finite  Element 
Method  for  Heterogeneous 
Microstructures 


Summary 

The  Voronoi  Cell  finite  element  model  (VCFEM)  has  been  successfully  developed  for 
materials  with  arbitrary  microstructural  distribution.  In  this  method,  the  finite  element 
mesh  evolves  naturally  by  Dirichlet  Tessellation  of  the  microstructure.  Composite  VCFEM 
for  small  deformation  plasticity  has  been  developed  by  expressing  the  element  stresses  in 
terms  of  polynomial  expansions  of  location  coordinates.  Though  this  works  well  for  discrete 
composites  with  inclusions,  its  effectiveness  diminishes  sharply  for  porous  materials  with 
voids.  The  effect  worsens  sharply  with  voids  of  arbitrary  shapes.  To  overcome  this  limitation, 
a  new  way  of  defining  stress  functions  is  introduced  in  this  work.  Based  on  a  transformation 
method  similar  to  the  Schwarz-Christoffel  conformal  mapping,  it  introduces  reciprocal  stress 
functions  that  are  derived  to  incorporate  shape  effects.  Several  numerical  experiments  are 
conducted  to  establish  the  strength  of  this  formulation.  The  effect  of  various  microstructural 
morphologies  on  the  overall  response  and  local  variables  are  studied. 


3.1  Introduction 

Increased  use  of  many  advanced  heterogeneous  materials  in  the  aerospace,  automotive  and 
electronics  industries,  has  prompted  widespread  research  for  understanding  their  deforma¬ 
tion  and  failure  mechanisms  in  the  recent  years.  Important  classes  of  materials  among  them 
are  metal/alloy  systems  with  microscopic  precipitates  and  pores,  and  composites  containing 
a  dispersion  of  fibers,  whiskers  or  particulates  in  the  matrix.  The  influence  of  microscopic 
heterogeneities  on  the  overall  behavior,  depends  on  morphological  characteristics  like  size, 
shape,  orientation  and  spatial  distribution  of  constiuent  phases,  as  well  as  on  their  proper¬ 
ties.  For  example,  ductility  is  reduced  with  increasing  volume  fraction  of  reinforcements  in 
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metal-matrix  composites.  Christman  et.  al.  [1]  have  shown  that  local  plastic  flow,  which 
yields  weaker  areas  of  inhomogeneous  deformation,  is  highly  sensitive  to  shape  of  the  rein¬ 
forcements  for  identical  global  stresses.  Studies  reflecting  the  details  of  actual  (not  idealized) 
microstructure,  are  indispensable  for  establishing  microstructure-property  relationships.  It  is 
through  these  studies,  that  the  effect  of  morphology  on  evolving  state  variables  like  stresses, 
plastic  strains,  void  growth,  and  material  variables  like  strain  hardening  and  flow  stress,  can 
be  established. 

A  number  of  micromechanical  theories  have  been  developed  to  predict  overall  constitu¬ 
tive  response  by  employing  continuum  mechanics  principles  at  the  microscopic  level.  Notable 
among  various  analytical  micromechanical  models  are  those  based  on  variational  approach 
using  extremum  principles  [2] ,  probabilistic  approaches  [3] ,  self-consistent  schemes  [4,  5]  and 
the  generalized  self  consistent  models  [6].  Though  analytical  micro-mechanical  models  are 
reasonably  effective  in  predicting  equivalent  material  properties  for  simple  geometries  and 
low  second  phase  volume  fractions,  arbitrary  distribution  of  shapes,  sizes  and  location  in 
real  materials,  cannot  be  deterministically  treated.  Constitutive  response  of  the  constituent 
phases  are  also  restricted,  and  predictions  with  large  property  mismatches  are  not  reliable. 
Despite  the  progress  in  analytical  modeling  of  brittle  heterogeneous  materials,  relatively  little 
has  been  done  for  ductile  materials.  Important  contributions  have  been  made  by  Tandon  and 
Weng  [7],  who  have  used  the  Mori-Tanaka  mean-stress  theory  in  conjuction  with  small  strain 
elastic-plastic  deformation  theory,  and  by  Dvorak  and  Bahei-El-Din  [8]  using  a  “vanishing 
fiber  diameter”  model,  for  rate-independent  plastic  matrix  with  elastic  fibers.  More  recently, 
Dvorak  and  Bahei-El-Din  [9]  have  developed  the  bimodal  plasticity  theory  for  nonhardening 
matrices.  Aboudi  [10]  has  predicted  the  response  of  composites  with  viscoplastic  matrices 
and  elastic  inclusions  and  voids.  Nemat-Nasser  and  coworkers  [11]  have  used  Fourier  series 
expansions  for  periodic  microstructures  to  develop  elastic-plastic  and  creep  models.  The 
appropriateness  of  these  models  for  many  real  materials  is  however  questionable.  This  is 
because  when  plastic  flow  occurs  in  the  ductile  phase,  the  deformation  is  no  longer  homo¬ 
geneous.  Local  properties  become  stress  dependent  and  the  overall  constitutive  response  is 
influenced  by  the  distributions  and  shapes  of  second  phase. 

Intractability  of  analytical  models  have  necessitated  the  introduction  of  Unit  Cell  for¬ 
mulations  [1,  13,  14]  using  computational  techniques  like  the  finite  element  method.  These 
models  generate  overall  material  response  through  detailed  discretization  of  a  representa¬ 
tive  material  element  (RME).  A  majority  of  these  models  make  assumptions  of  “perfect” 
local  periodicity  and  uniform  second  phase  distribution.  Effectively,  the  local  periodicity 
assumption  can  reduce  the  RME  to  a  basic  structural  element  (BSE),  thereby  making  the 
unit  cells  very  simple.  A  basic  structural  element  is  defined  as  the  smallest  element  of  the 
microstructure  reflecting  the  basic  geometrical  features  (see  figure  1(b)).  It  is  generally 
observed  though,  that  real  microstructures  rarely  possess  such  ordered  structures  as  are  nec¬ 
essary  for  this  idealization.  For  accurate  depiction  of  the  microstructural  evolution,  the  size 
of  conventional  FEM  elements  should  be  one  to  two  orders  of  magnitude  smaller  than  the 
BSE  size.  This  results  in  the  creation  of  very  large  meshes  and  enormous  computations, 
for  predictions  on  morphological  changes  and  their  effects  on  microstructural  failure.  This 
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limitation  of  conventional  finite  element  methods  necessitate  novel  computational  methods 
for  increased  effectiveness  in  micro-mechanics  studies. 

To  overcome  the  shortcomings  of  conventional  Unit  Cell  models  for  arbitrary  microstruc- 
tural  distribution,  Ghosh  and  coworkers  have  innovated  a  material  based  Voronoi  Cell  Finite 
Element  Method  (VCFEM)  [26,  27,  28,  29,  30,  19,  20].  In  this  method,  the  finite  element 
mesh  evolves  by  Dirichlet  Tessellation  of  the  microstructure  (see  [19],  resulting  in  a  net¬ 
work  of  multi-sided  “Voronoi”  polygons  or  cells  containing  one  heterogeneity  at  most  (figure 
3.1).  Effectively,  each  Voronoi  cell  element  may  be  identified  as  a  basic  structural  element 
of  the  microstructure.  In  VCFEM  formulations,  each  cell  with  the  embedded  heterogeneity 
is  treated  as  an  element  without  any  further  discretization.  It  can  thus  drastically  reduce 
preprocessing  efforts  as  well  as  actual  computations  expended  for  complex  RME’s.  Addition¬ 
ally,  it  has  been  demonstrated  by  Richmond  and  coworkers  [22]  tessellation  methods  can  be 
used  very  effectively  in  quantitative  characterization  of  micrographs.  Thus  VCFEM  can  be 
used  effectively  to  provide  a  direct  link  between  quantitative  metallography  and  mechanical 
response.  In  the  VCFEM  fomulations  for  small  deformation  plasticity  [26,  27],  the  stress 
field  within  each  element  is  expanded  as  polynomials  of  position  coordinates.  Though  this 
works  well  for  discrete  composites  with  inclusions,  its  effectiveness  diminishes  sharply  for 
porous  materials  with  voids.  Large  number  of  terms  are  required  in  the  stress  expressions 
to  yield  acceptable  accuracy.  The  accuracy  drops  sharply  with  voids  of  arbitrary  shapes. 
This  is  because  the  stress  functions  are  required  to  meet  the  zero  traction  condition  at 
the  void  boundary  and  undergo  very  large  gradients  near  the  interface.  For  composites,  this 
is  less  drastic  and  a  reasonable  number  of  terms  can  accommodate  the  traction  discontinuity. 

In  this  work,  a  new  method  of  defining  stress  functions  is  introduced  as  a  remedy.  Based 
on  a  transformation,  similar  to  the  Schwarz-Christoffel  conformal  mapping,  reciprocal  stress 
functions  are  derived  to  incorporate  shape  effects.  These  modifications  make  it  very  suitable 
for  application  to  a  wide  variety  of  heterogeneous  materials.  The  chapter  begins  with  a 
brief  overview  of  the  VCFEM  formulation  and  elucidates  the  limitations  of  pure  polynomial 
based  stress  expansions.  It  then  introduces  shape  based  stress  interpolations  and  solves 
several  numerical  experiments  to  establish  the  strength  of  this  formulation.  Finally,  the 
effect  of  various  microstructural  morphologies  on  overall  response  and  local  evolution  are 
qualitatively  studied. 

3.2  The  Voronoi  Cell  Finite  Element  Model  for  Heterogeneous 
Microstructures 

Voronoi  cells  resulting  from  Dirichlet  tessellation  of  a  heterogeneous  microstructure  make 
rather  unconventional  elements,  due  to  the  arbitrariness  in  the  number  of  edges.  The  appli¬ 
cation  of  conventional  displacement-based  finite  element  methods  to  these  elements,  suffer 
from  difficulties  associated  with  interelement  displacement  compatibility  and  rank  deficien¬ 
cies  of  the  stiffness  matrix.  The  Voronoi  Cell  finite  element  model  developed  by  Ghosh  and 
coworkers  [26,  27,  28,  29,  30]  avoids  these  difficulties  by  invoking  the  assumed  stress  hybrid 
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method  introduced  by  Pian  [33].  In  this  formulation,  independent  assumptions  are  made 
on  an  equilibriated  stress  field  in  the  interior  of  each  element  and  a  compatible  displace¬ 
ment  field  on  the  element  boundary.  Small  deformation,  elastic-plastic  analysis  of  materials, 
embedded  with  second  phase  inclusions,  has  shown  significant  promise  with  this  method 
[26,  27].  The  developments  use  an  assumed  stress  hybrid  formulation  for  rate  independent 
J2  flow  theory,  based  on  a  generalized  Hu-Washizu  principle,  originally  proposed  by  Atluri 
[34,  35]).  Though  the  details  of  this  method  are  available  in  [26,  27],  a  brief  overview  of  the 
formulation  essentials  is  presented  here  for  completeness. 

3.2.1  Variational  principle  with  assumed  stress  hybrid  method 

Consider  a  typical  representative  material  element  (RME)  11,  as  shown  in  figure  1.  The 
heterogeneous  domain  11  is  tessellated  into  N  Voronoi  cells  based  on  the  location,  shape 
and  size  of  N  heterogeneities  as  explained  in  [19]  .  The  matrix  phase  in  each  Voronoi 
cell  is  denoted  by  !!„  and  the  heterogeneity  (void  or  inclusion)is  denoted  by  flc.  Each 
element  boundary  dOe  is  assumed  to  be  comprised  of  three  mutually  disjoint  parts,  viz. 
(a)  prescribed  traction  boundary  Ft,,,,  (b)  prescribed  displacement  boundary  r„m,  and  (c) 
interelement  boundary  F^,  i.e.  dUe  =  Fi^UrtimUrm-  The  matrix-second  phase  interface 
dQc  has  an  outward  normal  n'",  while  is  the  outward  normal  to  dfig.  An  incremental  finite 
element  formulation  is  invoked  to  account  for  the  evolutionary  constitutive  equations  in  rate 
independent  plasticity.  At  the  beginning  of  the  p-th  increment,  let  <r  be  an  equilibriated 
stress  field  with  a  strain  field  e{<T,load  history),  and  u  be  a  compatible  displacement  field 
on  the  element  boundary.  Also  let  A<r  correspond  to  an  equilibriated  stress  increment  in  fig, 
Au  to  a  compatible  displacement  increment  on  dftei  and  At  to  a  traction  increment  on  Ttm- 
The  incremental  problem  is  solved  by  using  a  two  field  assumed  stress  hybrid  variational 
principle,  derived  from  an  element  energy  functional  as: 

ne(Ao-,  Au)  —  -  [  AB(<t,A<t)  dQ  -  f  e  :  Aa  dfl 
+  f  (o-  -t-  Ao-)  ■  n"  •  (u  -h  Au)  dCl-  f  (t  +  At)  •  (u  -f  Au)  dT 

JTtm 

-I  +  A<t^  -  -  A^^)  ■  ■  {u' +  Au')  dCl  (3.1) 

where  Au'  the  displacement  of  the  interface  and  AB  is  the  increment  in  element  compli¬ 
mentary  energy.  Superscripts  m  and  c  represent  respectively  the  matrix  and  second  phase 
parts  of  the  Voronoi  cell  element.  The  energy  functional  for  the  entire  domain  is  obtained 
by  adding  each  element  functional  as 


n  =  Ene  (3.2) 

e=l 

The  first  variation  of  Ffg  with  respect  to  the  stress  increments  A(r,  results  in  the  kinematic 
relations  as  the  Euler  equation. 


VAu  =  Ae  in  fig 


(3.3) 
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Figure  3.1:  (a)  Representative  material  element  (b)  Voronoi  ceU  finite  element 

while  the  first  variation  of  11  with  respect  to  boundary  displacement  increments  Au  yields 
traction  conditions  as  Euler  equations, 

(<7  +  A<t)  •  n®'''  =  — (<T  +  A<t)  •  n®“  on  Interelement  traction  reciprocity 

{<T  +  A<t)  •  n®  =  t  +  At  on  Ft^  Traction  boundary  conditions 

(cr'^  +  A<r‘')  •  =  (rr”*  +  Ao-”*)  •  n‘^  on  dVlc  Interface  traction  reciprocity  (3.4) 

Equilibriated  stress  increments  A<r,  constitutive  relations,  along  with  the  incremented  form 
of  the  energy  functional  completely  define  the  p-th  increment  problem. 


3.2.2  Element  formulations  and  assumptions 

In  the  Voronoi  cell  finite  element  model  (VCFEM)  formulation,  independent  assumptions  on 
the  equilibriated  incremental  stress  field  in  ffe  and  on  the  compatible  incremental  displace¬ 
ment  field  on  dVl  are  made.  Furthermore,  independent  assumptions  on  stress  increments  A«t 
are  made  in  the  matrix  and  heterogeneity  phases  to  accommodate  stress  jumps  across  the  in¬ 
terface.  Use  of  Airy’s  stress  functions  $(a;,p)  is  a  convenient  method  of  deriving  polynomial 
forms  of  stress  increments  in  two-dimensional  analysis,and  is  expressed  as: 

~  ~  ’  ^^^y  ~ 

Different  expressions  may  be  assumed  for  $  in  the  matrix  and  inclusion  phases.  For  example, 
a  third  order  complete  polynomial  expansion  of  may  be  assumed  to  give  the  stress 
increments  in  the  matrix  phase  as 

0  0  0  X  0 

1x00  y 
0  0  1  —y  —X 


A(7"  1 

1  J 

^  - 

1  y 
0  0 
0  0 


dxdy 


(3.5) 
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where  {A^”‘}’s  correspond  to  a  set  of  yet  undetermined  stress  coefficients.  In  general,  the 
stress  functions  can  be  arbitrary  functions  of  location,  yielding  stress  increments  in  the  form 


{  Ait”  } 


P”(t,!/)]{  A/3”  } 


{  } 


P'(i,!/)  ]{  A/3'  } 


(3.7) 


Compatible  displacement  increments  may  be  generated  by  interpolation  in  terms  of  gener¬ 
alized  nodal  values  on  the  element  boundary  dflg  as  well  as  on  the  interface  dilc-  Note  that 
compatible  displacement  fields  on  dQc  imply  perfect  bonding  in  the  case  of  composite  ma¬ 
terials.  A  linear  interpolation  of  displacement  increments  on  the  fth  and  (i  +  l)th  boundary 
nodal  points,  for  example  may  be  expressed  as, 


{Au} 


Aux 

AUy 


Aq2i-i 

o 

o 

! 

Aq2t 

0  1  —  s/k  0  s/li  ^ 

Aq2i+i 

.  Aq2i+2  j 

where  k  is  the  length  of  segment  and  s  is  the  distance  from  the  i-th  node.  The  displacement 
increments  on  the  element  boundary  and  interface  may  be  written  as. 


{Au}  =  [L^]{Aq} 


(Au'}  =  [L‘=]{Aq'}  (3.9) 

where  {Aq}  and  {Aq'}  are  generalized  displacement  increment  vectors.  Substituting  el¬ 
ement  approximations  for  stresses  (3.7),  and  displacements  (3.9),  in  the  energy  functional 
(3.1),  and  setting  the  first  variations  with  respect  to  the  stress  parameters  A/3™  and  A/3‘^ 
respectively  to  zero,  results  in  the  following  two  weak  forms  of  the  kinematic  relations  (3.3), 

/  [P™ne -f  Ae}df)  =  /  [P™]^[n^][L^]dn  {Aq}  -  [  [P™]^[n‘=][L^]df]  {Aq'} 

JdQe  JdQc 

[  [P‘f{e  +  Ae]dSl  =  f  (Py[n'l[L'I<i!l  {Aq')  (3.10) 

JQc  JdQc 


For  two-dimensional  problems  [n^]  and  [n^^]  are  (3  x  2)  matrices,  consisting  of  the  normal 
components  in  the  form: 


[n(x)] 


^^/(x)  0 

0  nj,(x) 
TZy(x)  772,(x) 


Setting  the  first  variation  of  the  total  energy  functional  (3.2)  with  respect  to  Aq  and 
Aq'  to  zero  results  in  the  weak  form  of  the  traction  reciprocity  conditions  as 


E 


e=l 


+  A/3™  \  _ 
f3^  +  A^^  j- 
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/r™lL*F{t  +  At}ifi 
{0} 


(3.11) 


N 

E 

e=l 


For  an  elastic-plastic  material,  the  strain  increments  Ae  in  equation  (3.10)  are  non-linear 
functions  of  the  current  state  of  stress  <r  as  well  as  of  their  increments  A<t.  The  non-linear 
finite  element  equations  (3.10)  and  (3.11)  are  solved  for  the  stress  parameters  (A/3”,A^'^) 
and  the  nodal  displacement  increments  (Aq,  Aq')  at  the  p-th  increment. 


3.2.3  Constitutive  Relations 

A  rate  independent  small  deformation  elastio-plastic  constitutive  relation,  following  J2  flow 
theory  with  isotropic  hardening,  is  considered  in  this  work.  A  brief  account  of  the  numerical 
integration  of  the  constitutive  relations  is  presented  here.  More  details  are  given  in  [27]  and  a 
similar  method  is  also  presented  in  [26].  An  additive  decomposition  of  the  strain  increments 
Ae  into  an  elastic  part  Ae®  and  a  plastic  part  Ae^^  is  assumed,  i.e. 

Ae  =  Ae®  -(-  Ae^'  (3.12) 

The  elastic  part  of  the  strain  increment  tensor  is  obtained  by  the  inner  product  of  the 
compliance  tensor  with  the  stress  increment  tensor.  The  yield  surface  in  stress-space  at  the 
beginning  of  the  p— th  increment  is  expressed  as, 

FP(€,Ae)  =  y|o^o-'  (3.13) 

where  <t'  is  the  stress  deviator  and  Y (e,  Ae)  is  the  radius  of  the  flow  surface.  The  plastic 
strain  increment  is  obtained  by  numerically  integrating  the  flow  rule  by  the  backward 
Euler  method  to  yield  : 

AeP'  =  AA(or  +  Ao-)'  (3.14) 

where  AA  is  a  non-negative  incremental  flow  parameter.  Since  Ae  is  in  general  a  function 
of  AA,  <r  and  A<t,  the  flow  surface  radius  can  be  expressed  in  the  form, 

FP+i  =  yP+^(AA,o-,A<r) 

and  the  flow  parameter  AA  can  be  evaluated  from  the  following  relations. 

A  A  =  0  if  (elastic  unloading)  (3.15) 

if  >  Y^  (neutral  and  plastic  loading)  (3.16) 

Numerical  implementation  requires  computation  of  tangent  operators  through  linearized 
forms  of  the  constitutive  relations.  If  de  is  the  first  order  correction  to  the  current  strain 
increment  Ae,  and  dtr  is  the  corresponding  stress  increment  A<t  correction,  the  fourth  order 
elastic-plastic  compliance  tensor  (or  tangent  operator)  S  is  given  by  the  relation 

de  =  S  :  do¬ 


ll 


The  elastic  part  of  this  equation  expressed  as 


de^  —  Se  :  d(T 


(3.17) 


The  plastic  part  of  the  strain  correction  de^‘  however  requires  a  first  order  correction  dX  to 
the  current  flow  parameter  AA,  which  by  the  J2  flow  theory,  takes  the  form 


de^^  = 


9  (<T  +  A<t)' (8)  (o- +  Act)' 


:  da  =  S„/  :  da 


(3.18) 


where  H  is  &  linearized  hardening  modulus.  The  elasto-plastic  tangent  operator  S  is  obtained 
by  adding  equations  (3.17)  and  (3.18)  as 


S  —  Sj  +  Sj; 


(3.19) 


A  linearized  form  of  the  incremental  complementary  energy  functional  AB  in  equation  (3.1) 
can  now  be  expressed  in  terms  of  the  tangent  operator  as 


dB{a,da)  =  -da  :  S  :  da 


(3.20) 


Note  that  the  elastic-plastic  tangent  operator  S  in  equation  (3.19)  is  positive  definite  since  its 
components  are  individually  positive  definite.  Both  plane  stress  and  plane  strain  conditions 
have  been  solved  in  an  iterative  manner  by  this  basic  algorithms.  Details  are  provided  in 

[27] . 

3.2.4  Solution  Method  for  the  Weak  Form 

For  rate  independent  plasticity,  the  strain  increments  Ae(  A<t,  a)  are  nonlinear  functions  of 
the  stress  parameters  A/3™  and  A/3^  An  iterative  solution  process  is  therefore  invoked  for 
equation  (3.10)  to  evaluate  the  stresses,  given  the  nodal  displacement  increments  {Aq}  and 
{Aq'}.  Let  {d/3}“  correspond  to  the  correction  to  the  value  of  A/3  in  the  f-th  iteration,  i.e. 

{Ay3™}  =  {A/3™}‘ +  {(i/3™}‘ 

{A/3'=}  =  {A^y  +  {dfiy 

The  kinematic  equation  (3.10)  may  then  be  linearized  with  respect  to  A(3  to  yield; 


0  1  r  d/3^  y 


0  H. 


G.  -G, 

0  G^f 


q  +  Aq 
q'  +  Aq' 


lajP^ne  +  AeYdal 


In  a  condensed  form  this  can  be  restated  as: 


[n]{d^Y  =  [G]{Au}- 


/Q.[P"'r{^  +  Ae}'dD 
Aep'dQ 


V  e  =  1  •  •  •  iV 


(3.22) 
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where, 


[H™l  =  /  [p^i^isiip”].;^ ,  [Hj  =  /  [p'l^isiip'i-in 
[G,]  =  f  [P“l^[n*][L'l<in  ,  [G,„]  =  /  |P“|^[n'llL']in  ,  [G„]  =  /  [P']’'[n=l[L1rffi 

JdQe  JdQc  JdUc 

[S(x,  y)]  is  the  instantaneous  elastic-plastic  tangent  compliance  tensor  as  derived  in  equation 
(3.19).  A  quasi-Newton  iterative  solution  procedure  is  used  to  solve  equation  (3.22).  Details 
of  numerical  integration  schemes  for  plane  stress  and  plane  strain  conditions  are  presented 
in  [26,  27]. 

The  above  procedure  of  solving  for  the  stresses  take  place  within  an  iterative  loop,  in 
which  the  traction  reciprocity  conditions  (3.11)  are  solved  for  the  nodal  displacement  incre¬ 
ments  {Aq}  and  {Aq'}.  Proceeding  in  the  same  way  as  for  stresses,  let  {dq}-'  correspond 
to  the  correction  in  {Aq}  in  the  j—th.  iteration  of  (3.11),  i.e. 

(Aq)  =  {Aq}^  +  {dq}^' 

{Aq'}  =  {Aq'}^' +  {dq'}^ 

Substituting  equation  (3.22)  in  the  linearized  global  traction  reciprocity  equation  (3.11), 
with  respect  to  {Aq},  yields  the  matrix  equation: 

S[GnH]-‘[G]|  y =e{  |- 

^  [  /3n.[Vl’'[n']’-(P“l<in  0  1  /  ^”  +  A/3"  1' 

.^[-/9t..[Vf|n=f[P"]<JS!  V[L"|^M^(P']<in  J  \  ;3‘  +  A,3'  / 

or,  in  standard  finite  element  notation 

f;[K.niu}i  = 

e=l  e=l 


vf  /i>o.MVnP"l«  0  1  / /3"  +  A/3"  I  ’ 

.^[-VlVFMqP"]«  [P«|rffl  J  \  r  +  A/3'  / 


(3.23) 


With  known  traction  increments  on  Ttm  and  displacement  increments  on  the  linearized 
global  traction  reciprocity  condition  (3.23)  is  solved  iteratively  using  the  quasi-Newton 
method  for  nodal  displacement  increments. 


3.2.5  Numerical  Aspects  of  Voronoi  Cell  FEM 

Matrices  [Hm]  and  [He]  need  to  be  inverted  for  evaluating  the  element  stiffness  matrix  [Kg]. 
However,  Airy’s  stress  functions  represented  by  simple  polynomial  expressions  of  Cartesian 
coordinates  can  sometimes  lead  to  a  bad  conditioning  number  for  [H^]  and  [He].  This  can 
lead  to  considerable  numerical  inaccuracy  in  the  resulting  stiffness  matrix.  To  circumvent 
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this  problem,  Ghosh  et.al.  [26,  27]  have  used  scaled  polynomial  stress  functions  in  terms 
local  element  coordinates  (<f,  rj).  The  mapping  from  the  {x,  y)  coordinate  system  is  expressed 
as: 


7  =  {y-yc)li 

(3.24) 

where  (xc,  yc)  are  the  centroidal  coordinates  of  a  Voronoi  cell  element  and  /  is  a  scaling 
parameter  given  as: 


max{x 


Xc)  max{y-yc)) 


V  {x,y)  6 


The  scaling  to  ((f,  rj),  leads  to  an  approximate  range  of  variation  —  1  <  ^  <  1  and  —  1  <  <  1 

in  most  Voronoi  cell  elements.  Note  that  this  range  is  exactly  true  for  square  elements.  As 
an  example,  a  third  order  complete  polynomial  Airy’s  function  in  terms  of  (^,  rj)  gives  rise 
to: 


[P  = 


1  7/  0  0  0  ^  0 

0  0  1^00  y 

0  0  0  0  1  -77 


Another  factor  that  contributes  to  the  convergence  of  the  method  is  accurate  domain 
integration  to  evaluate  [H^]  and  [H^j.  For  polynomial  stress  functions,  the  matrices  are 
also  polynomial  functions  of  scaled  coordinates  (^,7/).  For  numerical  domain  integration  of 
these  functions,  the  Voronoi  cells  are  subdivided  into  quadrilaterals  (see  figure  3.4(a))  to  be 
mapped  into  square  master  domains.  Integration  is  then  performed  by  Gaussian  quadrature 
rule,  with  number  of  Gauss  points  determined  from  the  order  of  terms  in  [P]- 


3.3  Shape  Based  Stress  Representations 

An  important  criterion,  affecting  the  convergence  of  multiple  phase  Voronoi  cell  elements, 
is  the  proper  representation  of  stress  fields  in  each  of  the  constituent  phases.  In  general, 
choosing  matrix  stress  functions  from  micromechanics  considerations,  adds  considerably  to 
the  element  efficiency.  Three  different  conditions  that  are  indispensable  in  this  regard  are: 

1.  Stress  functions  should,  in  some  way,  account  for  the  shape  of  the  heterogeneity. 

2.  Effects  of  the  heterogeneity  shape  should  vanish  at  large  distances  from  the  interface, 
for  matrix  stress  functions  . 

3.  Shape  effects  in  matrix  stress  functions  should  facilitate  traction  reciprocity  at  the 
interface. 

For  both  composite  and  porous  materials,  the  first  two  considerations  imply  that  the  shape 
effect  should  be  dominant  near  the  interface,  but  vanish  in  the  far-held.  The  third  condition 
is  intended  to  counteract  interface  tractions  caused  by  the  inclusion  for  composites,  while 
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(b) 


Figure  3.2;  (a)  Actual  shape  of  a  Voronoi  Cell  with  embedded  heterogeneity  (b)  Transformed  shape 
of  the  VCE  by  approximate  Schwarz  Christoffel  transformation 


reduce  to  zero  interface  tractions  for  porous  materials. 

Pure  polynomial  forms  Airy’s  stress  functions  do  not  explicitly  account  for  the  shape  of 
the  heterogeneity,  and  consequently  requires  very  high  order  terms  (especially  for  porous 
media)  for  convergence.  The  following  example  in  figure  3.3(a)  is  considered  to  illustrate 
the  necessity  of  adequate  stress  field  representation.  A  single  unit  square  Voronoi  cell  el¬ 
ement,  with  a  circular  void  of  volume  fraction  Vf  =  2.5%  is  unidirectionally  stretched  by 
Auj;  =  0.008  units,  under  plane  stress  conditions.  The  material  is  elastic  with  Young’s 
modulus  E  =  250  and  Poisson  ratio  v  =  0.33.  Two  different  polynomial  stress  functions 
are  considered  for  solving  this  problem,  viz.  (a)  a  33  term,  7th  order  stress  function  ($) 
corresponding  to  a  5th  order  stress  representation  and  (b)  a  63  term,  10th  order  $  corre¬ 
sponding  to  a  8th  order  stress  representation.  The  deformed  configurations  are  depicted  in 
figure  3.3(b).  It  clearly  shows  that  lower  order  stress  representations  gives  rise  to  unsta¬ 
ble  deformation  modes  at  the  void  interface,  probably  trigerred  through  under-represented 
element  energy.  Extremely  poor  approximation  to  the  stress  concentration  is  observed  in 
figure  3.3(d),  which  shows  its  distribution  along  a  line  through  the  element  center.  Though 
the  interface  deformation  improves  with  higher  order  stress  terms  (figure  3.3(b)),  the  stress 
values  are  less  than  acceptable  near  the  interface.  The  results  emphasize  the  importance  of 
proper  matrix  stress  function  selection  in  VCFEM  formulation. 

Problems  of  stress  concentrations  around  voids  have  been  traditionally  handled  by  the  use 
of  specialized  stress  functions,  that  account  for  its  shape.  Analytical  solutions  by  Muskhe- 
lishvili  [27]  and  Sarin  [28]  use  the  Schwarz-Christoffel  conformal  mapping  to  transform  an 


15 


y/L 

(c)  (d) 


Figure  3.3:  (a)  VCFEM  mesh  for  a  unit  square  with  a  circular  hole  of  volume  fraction  Vf  =  2.5%  (b) 
Computed  deformed  shape  with  pure  polynomial  stress  functions  (P)  (c)  Computed  deformed  shape 
with  shape  based  stress  functions  (P+R)  (b)  Tensile  stress  distribution  along  the  line  x/L=0.5  for 
both  cases 


arbitrary  shaped  void  into  a  circle  in  the  complex  plane.  Analytic  functions  are  then  defined 
in  the  transformed  plane  and  used  to  generate  accurate  stress  functions.  Complex  functions 
and  conformal  mapping  techniques  have  been  used  by  Tong  et.al  [36]  and  Piltner  [37]  to  con¬ 
struct  trial  stress  functions  in  the  assumed  stress  hybrid  finite  element  method,  for  elastic 
problems  with  cracks  and  holes.  The  stress  functions  are  chosen  to  satisfy  the  biharmonic 
equation  because  of  linear  elastic  materials. 


Stress  functions  for  the  VCFEM  with  material  nonlinearity  are  based  on  the  three  con¬ 
ditions  mentioned  above.  Consider  a  a  typical  embedded  heterogeneity  as  shown  in  figure 
3.1(b).  Suppose  that  equation  of  the  interface  dQ,c  can  be  expressed  in  polar  coordinates  as 
^(r,  9)  =  0,  where  the  r  coordinate  is  measured  from  the  centroid  of  the  heterogeneity.  A 
Fourier  series  expansion  for  r  in  terms  of  the  polar  angle  6  may  be  expressed  as: 


r  =  Oo  +  ^  ctn  cos(n0)  +  ^  sin(n^)  on  dflc 

n  n 

where  the  Fourier  coefl&cients  an  and  6„  are  given  as: 

do  =  ^  /  rdO, 

27r  JdUc 

a„  =  —  /  rcos(nO)dCl  n  =  l,2,  ••• 

TT  JdQc 

bn  =  —  f  rsm(n6)dQ.  n  =  l,2,  ••• 

TT  JdQc 

The  interface  equation  may  then  be  expressed  from  (3.25)  as, 

gir,0)  =  f  -  —  -  XI  —  cos(n^)  -  X]  —  sin(n^)  =  0 

do  „  do  „  do 


(3.25) 


(3.26) 


Here  /  corresponds  to  a  function  that  transforms  any  arbitrary  shaped  interface  to  an  ap¬ 
proximate  unit  circle,  since  f{x,y)  =  1  on  dQ,c-  A  typical  numerical  transformation  of  a 
voided  cell  is  shown  in  figure  3.2.  The  mapped  function  f{r,9)  may  be  thought  of  as  a 
special  radial  coordinate  with  the  property  that  j  ^  0  as  (x,y)  oo.  This  function  is 
now  used  to  construct  shape  based  reciprocal  stress  functions  associated  with  polynomial 
functions.  Thus  if  the  polynomial  stress  function  is  represented  by  =  Y,p,q 

then  for  each  term  there  exists  reciprocal  terms  expressed  by  to  stresses 

equilibriating  the  traction  field  on  dQ,c-  This  is  written  as: 


_  Y'  _L  4-  .  .  .I 

!  \  Jp+q  ^  fp+q+1  ^  > 


P7<7 


The  Voronoi  cell  element  stress  function  in  the  matrix  phase  is  then  written  as: 


^7oiy 


+  $ 


m 

rec 


(3.27) 
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A4 

/\ty 


oo  1  oo  1  or 

P^Q  i=p-\~q  j  P,q  i=p+q  J 

oo  1  CO  1  /9f 

-Epr-'>)'(Aft,+  E  a/3„.-^)  +  e?'’>(«  E  P^s) 

P?9  ^=p+(?  J  P>9  i=p+g  / 


Coefficients  of  the  reciprocal  stress  function  A/3p,i,  in  first  set  of  terms  in  equation  (3.28), 
add  more  flexibility  to  polynomial  coefficients  A/?pg  for  matching  desired  tractions  at  the 
interface  {f  —  1).  The  gradient  of  /  in  second  set  of  terms  account  for  the  interface  shape. 
The  reciprocal  terms  in  equation  (3.28)  have  negligible  effects  on  the  traction  far  away  from 
the  heterogeneity,  by  /  becoming  extremely  large.  The  far-field  tractions  are  therefore  pro¬ 
duced  predominantly  by  polynomial  terms  in  the  stress  function  and  are  unaffected  by  the 
shape  of  the  heterogeneity.  Matrix  phase  stress  increments  are  obtained  from  the  traction 
increments  by  using  equation  (3.5)  as: 


=  [P^]{Ar}  (3.29) 

The  single  Voronoi  cell  element  problem  of  figure  3.3  is  solved  again  using  a  fourth  order 
complete  polynomial  expansion  for  the  stress  function  and  3  reciprocal  terms  per  polynomial 
term,  i.e.  p  +  q  —  2  •  •  •  4,  i  =  1  •  •  •  3.  This  corresponds  to  48  l3  parameters  in  the  matrix 
stress  field.  Figures  3.3(c)  and  (d)  shows  a  marked  improvement  in  the  deformed  void  shape 
as  well  as  in  the  stress  concentration  (P  +  R).  The  stress  field  enrichment  also  reduces  the 
number  of  required  /3  coefficients,  thus  enhancing  the  efficiency  significantly. 


3.3.1  Numerical  Implementation  of  Reciprocal  Stress  Functions 

In  practice,  smooth  heterogeneities  in  each  Voronoi  cell  element  are  modeled  as  discrete  n 
sided  polygons.  For  example,  figure  3.2(a)  shows  a  8-sided  polygonal  approximation  to  an 
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Figure  3.4:  Integration  subdomains  of  a  Voronoi  cell  for  (a)  pure  polynomial  stress  functions  (b) 
shape  based  stress  functions 


Figure  3.5:  Comparison  of  integration  accuracy  for  tensile  stress  distribution  with  old  and  new 
integration  domains 


elliptical  heterogeneity.  The  Fourier  coefficients  in  equation  (3.25)  are  generated  using  the 
values  of  {r,0)  at  interfacial  nodal  points.  Also  a  truncated  Fourier  series  is  used  for  the 
transformation  function  /  in  equation  (3.26).  Deviations  from  a  circle  in  figure  3.2(b)  cor¬ 
respond  to  these  numerical  aprroximations  in  the  Fourier  expansion. 

As  mentioned  in  section  2.5,  numerical  integration  of  [H„]  is  performed  by  Gauss  quadra¬ 
ture  rules,  applied  over  quadrilaterals  subdividing  Voronoi  cells.  For  polynomial  stress  func¬ 
tions,  this  can  be  done  with  accurately  by  estimating  the  desired  number  of  integration 
points  from  the  order  of  the  polynomial.  However,  when  reciprocal  stress  functions  are  used, 
the  number  of  quadrature  points  are  not  easily  defined.  Steep  stress  gradients  arising  from 
rapidly  changing  (i)  terms  occur  near  the  interface,  but  stabilize  in  the  far  field.  This  calls 
for  a  large  number  of  integration  points  for  quadrilaterals  in  figure  3.4(a).  Alternatively 
the  matrix  in  each  Voronoi  cell  is  divided  into  a  narrow  band  around  the  interface,  and 
the  remainder  of  the  cell,  as  shown  in  figure  3.4(b).  Both  these  regions  are  subdivided  into 
quadrilaterals  for  numerical  integration.  Quadrilaterals  in  the  band  are  integrated  using 
higher  number  of  integration  points  to  accommodate  the  sharp  gradients.  To  examine  the 
effectiveness  of  this  scheme,  the  voided  element  example  is  solved  again  using  original  and 
modified  integration  domains.  A  set  of  14x14  or  196  Gauss  points  per  quadrilateral  is  used 
for  the  old  scheme,  whereas  the  new  scheme  has  4x4  or  16  Gauss  points  for  the  inner  band 
and  3x3  or  9  Gauss  points  elsewhere.  The  resulting  stress  distribution  along  the  mid-section 
in  figure  3.5  indicates  that  the  modified  scheme  enjoys  similar  accuracy  for  a  considerably 
low  effort. 

3.3.2  A  Simple  Convergence  Study 

A  study  is  conducted  to  investigate  the  sensitivity  of  the  Voronoi  cell  elements  to  the  order 
terms  in  the  stress  function,  for  convergence.  Simulations  are  conducted  by  (a)  varying  the 
polynomial  order  with  constant  number  of  reciprocal  terms,  and  (b)  varying  the  order  of 
reciprocal  terms,  keeping  the  polynomial  terms  constant.  The  VCFEM  mesh  is  a  single 
element  with  a  heterogeneity  in  the  form  of  a  circular  voids  or  inclusion  as  shown  in  figure 
3.3.  It  is  subjected  to  uniaxial  stretching  with  Aui/T  =  0.008  under  plain  strain  conditions. 
Material  properties  are  as  follows: 

Matrix  (Elastic-Plastic) 

E  =  69  GPa  (Young’s  Modulus  ),  n  =  0.33  (Poisson’s  Ratio) 

Yq  =  Yi  MPa  (Initial  Yield  Stress),  a^qv  =  Fo  +  ^^eqv^  (Post  Yield  Hardening  Law  ) 
Inclusion  (Elastic) 

E  =  410  GPa  (Young’s  Modulus  ),  w  =  0.2  (Poisson’s  Ratio) 

Results  of  the  VCFEM  are  compared  with  those  generated  by  a  general  purpose  finite 
element  code  ABAQUS  with  an  extremely  refined  (approximately  1200  elements)  mesh  of 
QUAD4  elements.  The  first  example  is  for  a  circular  void  of  volume  fractions  Vf  =  5%  in  a 
unit  square  matrix.  3  different  ranges  for  the  matrix  stress  function  in  equation  (3.29)  are 
considered.  They  include: 

•  (a)  In  P  +  9  =  2. .4  ,  In  P  +  9  =  2. .4,  i  =  1..2  ;  36  terms  in  $ 
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Figure  3.6:  Tensile  stress  distribution  along  x/L  =  0.5  with  various  stress  functions  for  a  square 
domain  with  void,  at  void  volume  fraction  Vj  =  5% 

•  (b)  In  p-\-q  =  2..4  ,  In  p  +  9  =  2. .4,  i  -  1..3  ;  48  terms  in  $ 

•  (c)  In  p  +  q  =  2..6,  In  p  +  q  =  2..4,  i  =  1..3  ;  61  terms  in  $ 

Cases  (a)  and  (b)  are  intended  to  study  the  effect  of  additional  reciprocal  terms,  while  cases 
(b)  and  (c)  study  the  effect  of  the  polynomial  terms.  A  comparison  of  resulting  stresses  along 
the  mid-section  in  figure  3.6  shows  that  an  increase  in  the  number  of  reciprocal  terms  from 
case  (a)  to  case  (b)  significantly  improves  the  results.  On  the  other  hand,  an  increase  in  the 
number  of  polynomial  terms  from  case  (b)  to  case  (c)  does  not  significantly  alter  results. 
Thus,  from  both  efficiency  and  accuracy  point  of  view,  with  48  terms  (case  b)  appears 
to  be  the  optimal  choice  for  Voronoi  cells  with  voids. 

A  volume  fraction  of  Vj  =  35%  is  considered  for  the  convergence  test  of  a  composite 
cell  with  single  inclusion.  The  stress  field  in  the  inclusion  is  generated  with  a  25  term 
polynomial  stress  function  (i.e.  ^poiy-  p-f  9  =  2. .6)).  Different  cases  studies  for  the  matrix 
stress  function  include: 


•  (a)  In  p  +  9 

=  2..5  , 

In$™,:  p  +  9 

=  2,  i  =  1..3  ; 

24  terms  in  $ 

•  (b)  In$™,j,:  p  +  q 

=  2..6  , 

In$-,:  p-f9 

=  2,  i  =  1..3  ; 

34  terms  in  $ 

•  (c)  In  p  +  9 

=  2..6, 

In$™,:  p  +  9 

=  2..3,  i  =  1..3  ; 

46  terms  in  $ 

Cases  (a)  and  (b)  exhibit  the  effect  of  adding  polynomial  terms,  while  cases  (b)  and  (c)  study 
the  effect  of  reciprocal  terms.  Figure  3.6  shows  that  an  increase  in  the  number  of  polynomial 
terms  from  (a)  to  (b)  produces  a  significantly  improved  stress  distribution  in  the  inclusion. 
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domain  with  inclusion,  at  inclusion  volume  fraction  Vf  =  35% 


An  increase  in  the  number  of  reciprocal  terms  from  (b)  to  (c)  does  not  yield  any  significant 
advantages.  The  34  term  expansion  (case  b)  is  thus  chosen  as  the  optimal  stress  function 
for  Voronoi  cells  with  inclusions. 

3.4  VCFEM  for  Porous  Microstructures 

3.4.1  Numerical  Validation 

The  competence  of  the  Voronoi  cell  finite  element  model  in  analyzing  porous  microstructures 
is  ascertained  by  comparison  with  two  conventional  FEM  packages  (ABAQUS  and  ANSA'S). 
Assuming  that  the  domain  of  investigation  corresponds  to  a  microstructural  representative 
material  element  or  RME,  both  macroscopic  (volume  averaged)  and  microscopic  (detailed) 
behavior  are  investigated.  The  macroscopic  variables  are  denoted  with  an  overbar,  with 
macroscopic  stresses  &  expressed  as: 


<T  = 


(3.30) 


Regular  Packing 

Square  edge  packed  representative  material  elements  (RME)  with  circular  void  of  volume 
fractions  Vj  =  20%,  35%,  50%  and  square  shaped  void  of  volume  fraction  Vf  —  25%,  50%,  75% 
are  considered.  The  matrix  material  is  Aluminum  with  the  following  properties: 

E  -Q9  GPa  (Young’s  Modulus  ),  v  =  0.33  (Poisson’s  Ratio) 

Fo  =  43  MPa  (Initial  Yield  Stress),  aeqv  =  Fo  +  (Post  Yield  Hardening  Law  ) 

The  VCFEM  meshes  consisting  of  only  one  element  are  shown  in  figures  3.8(a)  and  (b),  while 
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the  ABAQUS  meshes  are  shown  in  figure  3.8(c)  and  (d).  Plane  strain  uniaxial  tension  loading 
is  considered  to  a  maximum  of  0.5%  macroscopic  strain.  Periodic  construct  of  the  RME  is 
achieved  through  repeatability  conditions  on  the  traction  free  face  {y  —  L).  This  constraint 
condition  requires  the  RME  to  deform  into  rectangular  shapes,  and  for  the  uniaxial  case  is 
written  as: 

Auy  =  Au  on  y  =  L  such  that  I  tydx  —  0 

Jy=L 

The  VCFEM  analysis  uses  a  48  term  stress  function  {p+q  =  2. .4  and  i  =  1..3)  in  the  matrix. 
The  VCFE  mesh  consists  of  8  nodes  on  the  element  boundary  and  8  nodes  on  the  interface. 
Thus  VCFEM  uses  a  total  80  D.O.F  (48  D.O.F  for  the  stress  parameters  and  32  D.O.F 
for  the  nodal  displacements,  in  contrast  to  the  ABAQUS  meshes  that  have  1972  and  1246 
D.O.F  for  the  circular  and  square  voids  respectively.  The  macroscopic  behavior  of  these  two 
RME’s  are  shown  in  figures  3.9(a)  and  (b),  and  shows  excellent  convergence  to  the  ABAQUS 
results.  Figures  3.10(a)  and  (b)  shows  the  stress  distributions  along  the  mid-section  of  the 
void  at  a  0.5%  macroscopic  strain  level.  These  results  establish  the  accuracy  of  the  VCFEM 
formulations  for  porous  materials  with  different  void  shapes. 

Random  Packing 

VCFEM  analysis  of  a  RME  consisting  of  29  randomly  located  circular  voids  (V/  =  20%) 
is  compared  with  results  generated  by  ANSYS.  Void  locations  are  generated  by  a  random 
number  generator.  Figures  3.11(b)  and  (a)  show  the  VCFEM  mesh  with  29  elements  and 
the  ANSYS  mesh  with  5282  QUAD4  elements  respectively.  Material  properties  and  loading 
conditions  are  identical  to  those  in  the  previous  example,  with  the  maximum  macroscopic 
strain  changed  to  0.8%.  Comparison  of  the  macroscopic  response  in  figure  3.12(a)  clearly 
establishes  VCFEM  as  an  accurate  method  for  modeling  overall  behavior  of  random  mi¬ 
crostructures.  The  microscopic  stress  distribution  through  the  section  A-A  in  figure  3.11(b) 
at  0.8%  strain  is  depicted  in  figure  3.12(b).  Once  again  the  comparison  is  very  satisfactory, 
with  VCFEM  producing  similar  patterns  and  peak  stresses  with  a  maximum  difference  of 
18.7%. 

The  equivalent  microscopic  plaatic  strain  distribution  at  a  macroscopic  strain  =  0.5%, 
is  shown  in  figure  3.13.  The  RME’s  include  square  edge  packing  with  circular  voids  (figure 
3.8(a))  with  different  volume  fractions  and  random  packing  (figure  3.11(b)).  As  noted  in  [31], 
figure  3.13(a)  shows  that  at  low  void  volume  fractions,  plastic  strains  are  localized  in  narrow 
diagonal  bands  in  the  form  of  ligaments  that  originate  from  the  void  boundary.  As  the  void 
boundary  approaches  the  cell  boundary  with  increasing  volume  fractions  (figures  3.13(a) 
and  (b)),  the  localization  effects  abate.  For  the  random  packed  RME  in  figure  3.13(d),  the 
maximum  plastic  strain  is  concenterated  around  a  few  void  boundaries.  It  is  interesting  to 
note  that  the  maximum  plastic  strain  for  the  random  RME  increases  approximately  tenfold 
with  respect  to  the  square  edge  packing  with  same  volume  fraction. 
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Figure  3.9:  Macroscopic  stress-strain  response  at  various  volume  fractions  for  square  edge  packed 
(a)  circular  voids  (b)  square  voids 


Microscopic  tensile  stress  (MPa)  Microscopic 


Microscopic  tensile  stress  (in  MPa) 


Figure  3.11:  (a)  ANSYS  and  (b)  VCFEM  meshes  for  random  packed  RME  with  circular  voids 
(Vf  =  20%) 


Figure  3.12;  (a)  Macroscopic  stress-strain  response  and  (b)  Microscopic  stress  distribution  at  0.8% 
strain  for  random  packed  RME  with  circular  voids,  along  xj L  ■=  0.5 
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3,4,2  Effect  Of  Microstructural  Morphology  On  Yield  Surfaces 

Effects  of  size,  distribution  and  shape  on  macroscopic  yield  surfaces  used  to  characterize 
the  onset  of  overall  elastic-plastic  behavior  are  considered  in  this  section.  The  Tvergaard- 
Gurson  (T-G)  model  [31]  for  square  edge  packed  cylindrical  voids,  is  used  as  a  benchmark  for 
comparison.  The  T-G  model  assumes  that  the  macroscopic  yield  surface  is  independently 
a  function  of  hydrostatic  and  equivalent  components  (first  and  second  invariants)  of  the 
macroscopic  stress  and  is  expressed  as: 

=  )^  +  29i  /cos/i(^^2— =  0 

0‘q  Z  (Jq 

where  aeqv  =  and  ahyd  =  ^kk,  ^ij  is  the  macroscopic  deviatoric  stress,  ao  cor¬ 

responds  to  the  matrix  flow  stress,  /  is  the  void  volume  fraction  and  coefficients  =  1 
and  q2  =  1.5.  Yield  surfaces  are  generated  from  VCFE  results  by  applying  various  differ¬ 
ent  boundary  conditions  on  the  RME  with  the  repeatibility  assumption.  These  conditions 
include  simple  shear  (tension  in  the  x-direction  and  compression  in  the  y-direction),  biaxial 
tension  (equal  tension  in  both  directions)  and  a  number  of  loading  types  in  between.  For 
each  of  these  loadings,  an  overall  equivalent  stress-strain  curve  is  drawn.  The  yield  point  is 
then  designated  as  the  point  where  the  tangent  modulus  undergoes  a  rapid  reduction.  This 
assumption  is  consistent  with  flnite  element  calculations  performed  by  Horn  and  McMeek- 
ing  [37]  for  spherical  voids.  For  this  problem,  the  matrix  material  is  assumed  to  have  the 
following  properties. 

E  =250  Units  (Young’s  Modulus  ),  u  =  0.33  (Poisson’s  Ratio) 

Yo  =  lUnit  (Initial  Yield  Stress),  =  0.00354  *  (cre^„)^°  (Post  Yield  Hardening  Law 


Effect  of  Void  Location 

For  circular  voids,  the  effects  of  three  position  distributions,  viz.  (a)  Square  edge  packing 
(S.E)  (figure  3.14(a))  (b)  Square  diagonal  packing  (S.D)  (figure  3.14(b))  and  (c)  Random 
packing  (figure  3.14(c)),  on  the  yield  surface  is  studied.  Four  different  void  volume  fractions 
of  2.5%,  5%,  7.5%  and  10%  are  considered  to  examine  the  dependence  on  volume  fraction. 
Figure  3.15  shows  the  yield  function  plotted  in  terms  of  equivalent  stress  as  a  function  of  hy¬ 
drostatic  stress  for  these  cases.  As  is  expected,  VCFEM  simulated  yield  functions  for  square 
edge  packing  are  in  very  close  agreement  with  the  Tvergaard-Gurson  predications  at  all  vol¬ 
ume  fractions.  It  is  interesting  to  note  that  the  yield  surfaces  for  random  packing  also  closely 
follow  the  Tvergaard-Gurson  model  at  lower  volume  fractions,  though  the  deviations  slightly 
increase  at  higher  volume  fractions.  The  square  diagonal  packing,  however,  is  at  significant 
variance  with  the  Tvergaard-Gurson  model  at  all  volume  fractions.  From  these  simulations, 
it  appears  that  hydrostatic  stress  has  a  less  pronounced  effect  on  the  yield  surface  for  the 
square  diagonal  packing.  The  effect  of  void  distribution  is  much  more  pronounced  at  higher 
volume  fractions  as  evidenced  from  figure  3.17,  where  yield  surfaces  are  depicted  at  Vy  =  35%. 


Size  Effect 

The  yield  surface  of  circular  voids  in  a  square  edge  packing  arrangements  at  volume  fractions 
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(c)  (d) 

Figure  3.14:  VCFEM  meshes  for  (a)  square  edge  packing  (SE),  (b)  square  diagonal  packing  (SD), 
(c)random  packing  (R)  ,  and  (d)  hexagonal  close  packing  (HCP)  RMEs 
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Figure  3.16:  Size  effect  on  yield  surfaces  of  square  edge  packed  RMEs’  with  circular  voids 


Figure  3.17:  Initial  macroscopic  yield  surfaces  at  a  higher  void  volume  fraction  (V)  =  35%) 


Figure  3.18:  VCFEM  meshes  for  elliptic  voids  in  (a)  square  edge,  and  (b)  random  packing  arrange¬ 
ments 

of  Vf  =  25%,  35%,  50%  are  shown  in  figure  3.16.  The  equivalent  stress  -  hydrostatic  stress 
plot  shows  that  the  yield  surface  shrinks  rapidly  with  increase  in  void  size.  A  comparison 
with  the  Tvergaard-Gurson  model  a.t  Vf  —  25%  shows  significant  variance  with  the  VCFEM 
results.  This  may  be  attributed  to  the  fact  that  the  continuum  derivation  of  Tvergaard- 
Gurson  model  is  only  valid  for  smaller  volume  fractions,  at  which  the  voids  do  not  interfere 
with  tractions  boundaries  of  the  RME. 

Void  Shape  and  Orientation  Effects 

Orientation  and  void  shape  effects  on  the  yield  surface  are  investigated  in  this  example. 
Figure  3.18(a)  shows  an  elliptic  void  in  a  square  edge  arrangement.  The  yield  surfaces  in 
figure  3.21  are  drawn  for  this  RME  with  a/b  =  1.5  with  the  major  axis  oriented  at  0°,  30'’, 
and  45°  to  the  x-axis.  The  yield  surface  is  symmetric  about  the  =  0)  line  for  the  45° 
orientation  due  to  reflective  properties.  However  for  the  0°  and  30°  orientations,  the  yield 
surface  exhibits  nonsymmetry  as  the  stress  state  progresses  from  hydrostatic  -  0)  to 

pure  shear  {^kk  =  0).  This  indicates  that  if  a  biaxial  tension  loading  is  replaced  by  a  biaxial 
compression  loading  of  equal  magnitude,  the  yield  points  will  not  reciprocate,  even  though 
the  overall  equivalent  stress  remains  unchanged.  It  implies  that  for  each  point  {akk,^eqv) 
on  the  yield  surface,  the  point  {-akk^^eqv)  may  not  belongs  on  the  yield  surface.  This  leads 
to  the  conclusion  that  a  representation  of  yield  surface  in  the  d^kk'^eqv  plane  is  no  longer 
appropriate  for  elliptical  voids  with  orientation.  One  remedy  is  to  represent  the  yield  surface 
in  the  space  of  invariants,  i.e.  akk-^eqv-  \^\  space.  Alternatively,  the  yield  surface  can  be 
represented  in  the  a^x-^yy  plane  as  shown  in  figure  3.19.  The  yield  surface  for  the  0°  and 
30°  orienatations  deviate  from  that  for  the  45°  orienatation.  Anisotropy  of  the  0°  and  30° 
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Figure  3.19:  Orientation  effect  on  yield  surfaces  for  square  edge  packed  elliptic  voids  in  cFxx-^: 
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Figure  3.20:  Void  Shape  effect  on  the  yield  surface  for  elliptic  voids  in  a  square  edge  packing  at  (a) 
Vf  =  2.5%,  and  (b)  Vf  =  10% 
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orientation  yield  surfaces  is  observed  from  their  differences  with  the  45°  yield  surfaces  in  each 
of  the  4  quadrants.  Figure  3.20(a)  and  (b)  illustrates  the  effect  of  void  shape,  by  considering 
three  shapes  viz.  a/6  =  1.0,  ajh  =  1.5  and  ajh  =  2.0  for  volume  fractions  Vj  =  2.5%  and 
V/  =  10%.  The  orientation  is  set  at  45°  for  symmetry.  The  increase  in  aspect  ratio  the 
results  in  a  lowering  of  the  yield  surface,  with  the  shape  effect  becoming  more  pronounced 
at  higher  volume  fractions.  As  a  final  example,  a  comparison  between  yield  surfaces  for  a 
perfectly  random  (distribution,  shape  and  orientation)  RME  with  39  voids  in  figure  3.18(b), 
and  a  RME  with  randomly  distributed  circular  voids  in  figure  3.11(b),  is  made.  Figure 
3.22  shows  the  yield  surfaces  for  Vj  —  2.5%  and  Vj  =  10%.  In  general,  the  circular  voids 
result  in  higher  values  of  yield  functions  at  lower  hydrostatic  stress  values.  The  difference 
is  considerably  magnified  at  higher  volume  fractions,  with  the  conclusion  that  size  plays  a 
more  important  role  than  shape  for  this  behavior. 


3.5  VCFEM  for  Composite  Microstructures 

A  number  of  examples  for  various  types  of  composite  microstructures  using  polynomial  stress 
functions  have  been  discussed  in  Ghosh  and  Moorthy  [27].  In  this  section  a  few  examples  are 
considered  to  demonstrate  the  advantage  of  using  reciprocal  stress  functions  in  the  Voronoi 
cell  element  formulation  with  embedded  inclusions. 

RME  with  Randomly  Distributed  Inclusions 

The  representative  material  element  illustrated  in  figure  3.11,  is  considered  for  an  Aluminum 
matrix  embedded  with  brittle  Boron  inclusions.  The  material  properties  are  as  follows: 
Aluminum 
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Figure  3.22:  Effect  of  void-shape  on  the  yield  surface  for  random  packed  RMEs’  at  Vj  =  2.5%  and 
Vf  =  10% 


E  —  Q9  GPa  (Young’s  Modulus  ),  v  —  0.33  (Poisson’s  Ratio) 

Fo  =  43  MPa  (Initial  Yield  Stress),  =  Yq  +  (Post  Yield  Hardening  Law  ) 

Boron 

E  =  420  GP A  (Young’s  Modulus  ),  u  ~  0.2  (Poisson’s  Ratio) 

The  RME  is  uniaxially  loaded  under  plane  strain  conditions  to  a  macroscopic  strain  of  0.5%. 
A  .34  term  matrix  stress  function  is  used  for  the  modified  formulation,  and  compared  with  a 
25  term  (p  +  ^  =  2  •  •  •  6)  polynomial  stress  function  in  the  old  formulation.  The  macroscopic 
response  of  the  RME  is  shown  in  figure  3.23(a),  while  figure  3.23(b)  depicts  the  microscopic 
stress  distribution  along  x  =  9.1  L  at  0.5%  overall  strain.  VCFEM  results  are  comapred  with 
ANSYS  using  a  very  refined  mesh.  VCFEM  simulations  with  both  polynomial  stress  func¬ 
tions  (P)  and  polynomial  with  reciprocal  stress  functions  (P+R)  are  in  excellent  agreement 
with  ANSYS  results  for  the  macroscopic  response.  However  the  microscopic  stress  distri¬ 
butions  in  figure  3.23(b)  establishes  the  superiority  of  the  modified  VCFEM  formulation 
beyond  doubt  for  composites.  The  stress  shows  very  good  convergence  to  ANSYS  results, 
with  a  maximum  difference  of  9.8%  in  the  peak  stress. 

Size  Effects  on  Composite  Yield  Stress 

Effect  of  various  fiber  arrangements  on  overall  strength  of  metal  matrix  composites  are 
investigated  by  VCFEM  and  compared  with  results  by  Zahl  et.  al.  [34].  Zahl  et.  al. 
[34]  have  used  ABAQUS  to  obtain  their  results  for  MMC’s  with  ductile  matrices  and  rigid 
inclusions.  The  parameter  chosen  to  be  examined  is  the  asymptotic  yield  stress,  defined  as 
cT.v  in  the  power  law  hardening  model  a  =  (T^(e/eo)’”.  Bao  et.al.[33]  have  concluded,  that 
for  fully  developed  plastic  zones  with  rigid  inclusions,  the  overall  post-yield  behavior  of  the 
composite  can  be  obtained  by  replacing  amat  with  a  different  in  the  matrix  hardening 
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law  a  =  CTynati^l ^o)'^ ■  The  matrix  material  in  this  example  has  the  following  properties. 

E  —  88.867  Units  (Young’s  Modulus  ),  i/  =  0.33  (Poisson’s  Ratio) 

(^mat  =  1  Unit  (Initial  Yield  Stress),  =  (7„a«(ee///0.001)°-^  (Post  Yield  Hardening 

Law  ) 

A  square  edge  (S.E),  square  diagonal  (S.D),  hexagonal  close  packing  (H.C.P)  and  a  random 
packed  RME,  as  shownin  figure  3.14,  are  uniaxially  loaded  with  repeatability  conditions. 
Figure  3.24  compares  the  VCFEM  generated  asymptotic  yield  stresses  with  results  in  [34], 
for  different  volume  fractions.  Excellent  agreement  is  obtained  for  all  packings. 

Comparisons  with  Some  Analytical  and  Experimental  Results 

Plane  strain  VCFEM  results  are  compared  with  analytical  2  and  3-phase  model  predictions 
due  to  Zhao  and  Weng  [35]  and  experimental  results  due  to  Adams  [42]  for  unidirectional 
composites.  The  material  properties  are: 

Aluminum 

E  =  58  GPa  (Young’s  Modulus  ),  u  =  0.33  (Poisson’s  Ratio) 

Yo  =  89  MPa  (Initial  Yield  Stress),  =  Yo  + 175(6^' MPa  (Post  Yield  Hardening 
Law  ) 

Boron 

E  =  385  GPA  (Young’s  Modulus  ),  i/  =  0.2  (Poisson’s  Ratio) 

A  2  x  1  rectangular  edge  packed  RME  with  a  circular  inclusion  of  V/  =  34%  is  loaded  in 
uniaxial  tension  with  repeatibility,  as  shown  in  figure  4.13(a).  Figure  4.13(b)  shows  that 
VCFEM  results  provide  a  better  agreement  with  experimental  results  than  the  analytical 
model.  This  is  probably  because  of  the  square  arrays  used  in  [35]  to  approximate  the  rect¬ 
angular  RME. 
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(a)  (b) 

Figure  3.25:  (a)  Rectangular  edge  packing  composite  RME  (V/  =  34%)  (b)Macroscopic  stress-strain 
responses  by  VCFEM,  analytical  and  experimental  methods 


3.6  Conclusions 

In  this  work,  the  Voronoi  cell  finite  element  model  is  further  developed  and  enhanced  for  ac¬ 
curate  and  efficient  2-D  analysis  of  heterogeneous  materials.  Specifically,  small  deformation 
elasto-plcistic  deformation  of  porous  and  composite  media,  with  an  arbitrary  dispersion  of 
microscopic  phases  is  considered.  Stress  functions  in  the  modified  scheme  are  motivated  from 
the  essential  characteristics  of  micromechanics  by  accounting  for  the  heterogeneity  shape  and 
influence.  This  is  achieved  by  transforming  arbitrary  shaped  heterogeneities  to  a  unit  circle 
using  a  method  similar  to  the  Schwarz-Christoffel  conformal  mapping.  Reciprocal  stress 
functions  of  singular  nature  are  then  constructed.  This  is  a  significant  advancement  over 
polynomial  stress  functions  commonly  used  in  hybrid  finite  element  methods.  It  enhances 
the  accuracy  for  various  microstructures  tremendously  at  very  moderate  computational  ef¬ 
forts.  It  would  take  a  prohibitively  large  number  of  terms  and  consequently  incur  enormous 
computational  costs  to  generate  results  of  comparable  accuracy  by  pure  polynomial  stress 
expansion. 

The  accuracy  and  efficiency  of  VCFEM  are  established  by  comparing  with  conventional 
FEM  commercial  packages.  For  a  wide  range  of  problems  VCFEM  delivers  very  similar 
accuracy  at  a  considerably  low  computational  effort.  This  is  evidenced  by  the  drastically 
reduced  degrees  of  freedom  needed  for  convergence,  compared  with  the  conventional  codes. 
The  D.O.F.  ratio  varies  from  as  low  cis  20  for  simple  microstructures,  to  even  100  for  more 
complex  cases.  This  translates  into  a  reduction  factor  of  15-30  in  the  CPU  time  for  execu¬ 
tion,  even  with  a  non-optimized  research  code.  Furthermore,  user  effort  required  to  generate 
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the  model  is  much  lower  for  VCFEM  than  for  commercial  packages. 

Mechanical  responses  and  their  qualitative  comparison  for  various  microstructures  are 
possible  through  VCFEM  analysis.  For  porous  materials,  VCFEM  generated  yield  functions 
are  compared  with  predictions  of  the  continuum  Tvergaard-Gurson  (T-G)  models.  Though 
VCFEM  results  agree  with  the  T-G  models  at  lower  void  volume  fractions,  it  is  noted  that 
increasing  void  size  produces  considerable  deviations.  Distribution  and  void  shape  also  have 
a  significant  influence  on  the  overall  performance,  and  these  effects  are  magnified  with  in¬ 
creasing  size.  This  is  also  evidenced  through  the  difference  in  yield  surfaces  generated  by 
VCFEM  and  the  T-G  model.  An  interesting  observation  is  that  void  orientation  plays  a  more 
important  role  than  void  shape  in  yielding  of  the  overall  material.  VCFEM  is  also  proved 
to  be  accurate  in  modeling  composite  microstructures.  VCFEM  comparisons  with  ANSYS 
and  ABAQUS  simulations,  as  well  as  with  experimental  and  analytical  results  establish  this 
claim. 

In  conclusion,  the  Voronoi  Cell  finite  element  model  has  emerged  as  an  important  tool 
for  analyzing  arbitrary  microstructures  in  many  materials.  It  strength  is  derived  from  the 
combination  of  basic  micromechanics  concepts  with  the  essentials  of  finite  element  methods. 
It  has  the  potential  for  direct  corelation  between  quantitative  characterization  and  material 
response.  Multiple  scale  coupling  between  structural  and  microstructural  scales  involving 
commercial  packages  at  the  structural  scale  is  easily  accomplished  with  this  method,  as 
shown  by  the  authors  in  [30,  21]. 
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Chapter  4 

Multiple  Scale  Modeling  of 
Heterogeneous  Materials  Using 
Asymptotic  Homogenization  and  the 
Voronoi  Cell  FEM 


Summary 

In  this  work,  a  multiple  scale  finite  element  model  (VCFEM-HOMO)  has  been  developed 
for  elastic-plastic  analysis  of  heterogeneous  (porous  and  composite)  materials  by  combining 
asymptotic  homogenization  theory  with  the  Voronoi  Cell  finite  element  model  (VCFEM). 
VCEEM  for  microstructural  modeling  originates  from  Dirichlet  tessellation  of  representative 
material  elements  at  sampling  points  in  the  structure.  Structural  modeling  is  done  by  the 
general  purpose  finite  element  code  ABAQUS,  and  interfacing  with  the  microscale  VCFEM 
analysis  is  done  through  the  user  subroutine  in  ABAQUS  for  material  constitutive  relation, 
UMAX.  Asymptotic  homogenization  in  UMAX  generates  macroscopic  material  parameters 
for  ABAQUS.  Following  the  macroscopic  analysis,  a  local  VCFEM  analysis  is  invoked  to 
depict  the  true  evolution  of  microstructural  state  variables.  Various  numerical  examples  are 
executed  for  validating  the  effectiveness  of  VCFEM-HOMO,  and  the  effect  of  size,  shape  and 
distribution  of  heterogeneities  on  local  and  global  response  is  examined. 


4.1  Introduction 

Xhe  last  three  decades  have  experienced  a  surge  in  the  advancement  of  science  and  technology 
for  heterogeneous  materials  that  contain  dispersions  of  multiple  phases  in  the  microstruc¬ 
ture.  Examples  are  metal/alloy  systems  with  second  phase  in  the  form  of  precipitates  and 
pores,  and  composite  materials  containing  a  dispersion  of  fibers,  whiskers  or  particulates  in 
the  matrix.  In  alloy  systems,  precipitates  and  pores  are  results  of  conventional  processing 
routes.  In  reinforced  composites,  stiff,  strong  and  brittle  second  phase  inclusions  e.g.  silicon 


46 


carbide,  boron  or  aluminum  oxide,  etc.  are  added  to  titanium,  nickel  or  aluminum  matrices 
to  enhance  flow  strength,  creep  resistance  and  wear  resistance  of  structures.  These  func¬ 
tionally  superior  materials  have  found  increasing  utilization  in  the  aerospace,  automotive, 
materials  and  ordnance  industries  for  replacing  some  of  the  traditionally  used  structural  ma¬ 
terials.  The  degree  of  mechanical  and  thermal  property  enhancements  depends  on  the  size, 
shape  and  properties  of  the  second  phase,  as  well  as  on  their  spatial  distribution  within  the 
matrix.  For  example,  ductility  is  reduced  with  increasing  volume  fraction  of  reinforcements 
in  metal-matrix  composites.  Christman  et.  al.  [1]  have  shown  that  local  plastic  flow  is 
highly  sensitive  to  shape  of  the  reinforcements  under  identical  global  stresses.  It  is  clear 
that  rigorous  fundamental  studies  are  warranted  for  understanding  deformation  and  failure 
mechanisms  prior  to  design  of  advanced  materials  in  high  performance  applications.  Stud¬ 
ies  reflecting  the  details  of  actual  heterogeneous  microstructures  are  indispensable  in  this 
respect.  It  is  through  these  models,  that  the  effect  of  second  phase  shapes,  sizes  and  distri¬ 
butions  on  evolving  state  variables  like  stresses,  plastic  strains,  void  initiation  and  growth, 
and  evolving  material  variables  like  strain  hardening  and  flow  stress  can  be  investigated. 

A  number  of  analytical  micromechanical  models  for  predicting  macroscopic  response  have 
evolved  within  the  framework  of  small  deformation  linear  elasticity  theory.  Notable  among 
them  are  models  based  on  :  (i)  variational  approach  using  extremum  principles  [2,  3],  (ii) 
probabilistic  approach  [4],  (iii)  self-consistent  schemes  [5,  6]  and  (iv)  the  generalized  self 
consistent  model  [7,  3].  These  models  follow  the  idea  of  equivalent  inclusion  methods  based 
on  eigenstrain  formulation.  Though  analytical  micro-mechanical  models  are  reasonably  ef¬ 
fective  in  predicting  equivalent  material  properties  for  relatively  simple  geometries  and  low 
volume  fractions,  they  are  often  incapable  in  depicting  the  evolution  of  stresses  and  strains 
in  the  microstructure.  Arbitrary  microstructural  morphology,  that  are  frequently  encoun¬ 
tered  in  actual  materials  cannot  be  deterministically  treated  with  these  models.  Constitutive 
response  of  the  constituent  phases  are  also  somewhat  restricted  and  predictions  with  large 
property  mismatches  are  not  very  reliable.  Additionally,  due  to  lack  of  proper  representation 
of  microscopic  stresses  and  strains,  they  cannot  capture  the  effect  of  local  inhomogeneities 
on  damage  and  failure.  The  state  of  the  art  in  analytical  modeling  for  ductile  materials  is 
not  as  mature.  Important  contributions  have  been  made  by  Tandon  and  Weng  [8]  for  small 
strain,  deformation  theory  of  elasto-plasticity  of  metals  with  elastic  particles,  and  by  Dvorak 
and  Bahei-El-Din  [9,  10]  for  rate-independent  plastic  matrix  and  elastic  fibers  using  a  “van¬ 
ishing  fiber  diameter”  model,  and  more  recently  the  bimodal  plasticity  theory.  Paley  and 
Aboudi  [11]  have  developed  a  semi-analytical  generalized  method  of  cells  for  elastic-plastic 
and  viscoplastic  materials,  while  Nemat-Nasser  and  coworkers  [12]  have  used  Fourier  series 
expansions  to  develop  elastic-plastic  and  creep  models.  The  applications  of  these  nonlinear 
models  to  complex  microstructures  in  many  real  materials  are  even  less  effective  than  the 
linear  models.  This  is  substantiated  by  the  fact  that  when  plastic  flow  occurs,  the  deforma¬ 
tion  is  no  longer  homogeneous.  Local  properties  become  stress  dependent  and  the  overall 
constitutive  response  is  influenced  by  the  distributions  and  shapes  of  second  phase. 

Intractability  of  analytical  models  in  situations  of  arbitrary  phase  dispersions  have  neces¬ 
sitated  the  introduction  of  Unit  Cell  formulations  [1,  13,  14]  using  computational  methods. 
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These  models  generate  overall  material  response  through  detailed  discretization  of  a  repre¬ 
sentative  volume  element  (RVE)  of  the  composite  microstructure.  In  unit  cell  models,  global 
properties  are  ascertained  by  assuming  macroscopic  periodicity  conditions  on  the  RVE’s. 
Despite  their  widespread  success,  the  unit  cell  models  are  faced  with  some  shortcomings  in 
modeling  real  materials,  stemming  from  global  and  local  periodicity  assumptions.  Global 
periodicity  assumes  that  the  same  RVE  is  valid  for  all  points  of  a  structure,  which  is  not 
realistic  for  many  commercial  materials.  A  majority  of  these  models  also  make  assumptions 
of  “perfect”  local  periodicity  and  uniform  distribution  of  the  second  phase.  Effectively,  this 
local  periodicity  assumption  can  reduce  the  representative  volume  element  (RVE)  to  a  basic 
structural  element  (BSE),  thereby  making  the  unit  cells  very  simple.  A  basic  structural 
element  is  defined  as  the  smallest  element  of  the  microstructure  reflecting  basic  geometrical 
features  (see  figure  1(c)).  Micrographs  of  many  real  materials  often  show  arbitrariness  in 
distribution  (figure  1(b))  thereby  making  this  assumption  too  restrictive.  Another  limitation 
is  that  unit  cell  models  typically  formulate  constitutive  relations  based  on  a  single  RVE  sub¬ 
jected  to  a  given  load  history.  This  makes  them  amenable  to  misrepresentation  of  the  actual 
interaction  between  macro-  and  micro-scale  deformations.  This  calls  for  simultaneous  com¬ 
putations  at  multiple  levels  for  concurrent  evolution  of  macroscopic  and  microscopic  variables. 

In  the  1970’s  a  mathematical  theory  called  Asymptotic  Homogenization  Theory  [15,  16], 
originated  for  analyzing  physical  systems  containing  two  or  more  length  scales.  It  is  suitable 
for  multi-phase  materials  in  which  the  natural  scales  are  the  microscopic  scale  characterized 
by  inter- heterogeneity/local  discontinuity  spacing  and  the  macroscopic  scale  characterizing 
overall  dimensions  of  the  structure.  This  method  is  based  on  asymptotic  expansions  of  dis¬ 
placement  and  stress  fields  about  their  respective  macroscopic  values  and  uses  variational 
principles  for  creating  a  link  between  scales.  Resulting  boundary  value  problems  at  the 
macroscopic  amd  microscopic  levels  are  solved  by  invoking  numerical  methods  at  both  lev¬ 
els.  Figure  1  shows  the  various  levels  in  a  structure  and  its  corresponding  computational 
model.  Computer  simulations  simultaneously  provide  global  response  through  homogenized 
material  parameters  and  averaged  stress/strain  fields,  and  microstructural  behavior  through 
a  depiction  of  local  stress/strain  fields.  This  method  can  overcome  shortcomings  inflicted 
by  global  periodicity  restrictions  in  unit  cell  models.  The  analysis  however,  makes  the  as¬ 
sumption  of  local  periodicity  through  the  introduction  of  spatially  repeated  microscopic  cells. 

The  finite  element  method  has  been  successfully  applied  in  conjunction  with  the  homog¬ 
enization  theory  for  analysis  of  linear  elastic,  particulate  and  fiber  reinforced  composites 
by  Toledano  and  Murakami  [17]  and  by  Guedes  and  Kikuchi  [18,  23].  These  models  are 
constructed  by  considering  very  simple  unit  cells  that  represent  uniform  distributions  in  the 
microstructure.  For  nonlinear  regimes,  homogenization  theories  have  been  established  by 
Suquet  [19].  Important  contributions  tothe  application  of  homogenization  method  in  finite 
deformation  elasto-plasticity  has  been  made  by  Guedes  [20]  and  Cheng  [21].  The  method 
has  been  extended  to  model  damage  evolution  by  fiber  rupture  for  linear  elastic  materials 
in  [25].  A  comprehensive  review  of  this  method  for  various  problems  is  presented  in  [22]. 
In  spite  of  its  promise,  the  homogenization  theory  has  not  enjoyed  a  wide  acceptance  thus 
far.  A  prime  reason  is  the  assumptions  of  over-simplified  microscopic  RVE’s  for  reducing 
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enormous  computing  costs. 


The  material  based  Voronoi  Cell  Finite  Element  Model  (VCFEM)  [26,  27,  28,  29,  30], 
developed  by  Ghosh  and  coworkers,  attempts  to  overcome  difficulties  in  modeling  arbitrary 
microstructures  by  conventional  finite  element  methods.  The  VCFEM  mesh  evolves  natu¬ 
rally  from  a  heterogeneous  microstructural  material  element  (RME)  by  Dirichlet  tessellation 
into  a  network  of  multi-sided  “Voronoi”  polygons.  The  polygons  contain  one  second  phase 
inclusion  each,  at  most  as  shown  in  figure  1(b).  A  robust  mesh  generator  to  create  these 
polygons  based  on  shape,  size  and  location  of  the  heterogeneities  is  developed  in  [30].  The 
multi-phase  Voronoi  polygons,  identified  with  the  basic  structural  elements  as  depicted  in 
figure  (Ic),  constitute  elements  in  VCFEM.  Element  formulations  have  been  developed  for 
micropolar  thermo-elasticity  problems  in  [29],  and  for  elastic-plastic  problems  in  [26,  27]. 
VCFEM  with  asymptotic  homogenization  for  elastic  problems  have  been  presented  in  [28]. 

This  work  presents  a  coupled  multiple  scale  computational  model  for  heterogeneous 
elastic-plastic  structures.  Only  two  dimensional  problems  are  considered.  Microstructural 
analysis  for  various  different  representative  material  element  arrangements,  is  done  with  the 
Voronoi  Cell  finite  element  model.  The  commercial  general  purpose  code  ABAQUS  [31]  is 
used  for  global  analysis  at  the  level  of  overall  structural  geometry  and  applied  loads.  The 
interfacing  between  macro-  and  micro-  calculations  are  done  through  the  user  subroutine 
window  UMAT  in  ABAQUS.  Numerical  examples  are  conducted  to  investigate  the  advan¬ 
tages  of  coupled  multiple  scale  analysis  over  other  unit  cell  and  continuum  based  theories. 
Effect  of  shapes,  sizes  and  locations  of  microscopic  inclusions  on  the  structural  performance 
and  the  evolution  of  microstructural  stresses  and  strains  are  also  studied. 


4.2  Asymptotic  Homogenization  for  Multiple  Scale  Analysis 

Consider  a  heterogeneous  body  occupying  a  region  il^tmcture  in  figure  4.1(a),  for  which  the 
microstructure  constitutes  of  spatially  periodic  representative  material  elements  (RME’s) 
as  shown  in  figure  4.1(b).  Only  small  deformation  elastic-plastic  deformation  of  the  body 
is  considered  in  the  absence  of  inertia.  The  body  is  subjected  to  a  system  of  time/load 
dependent  body  forces  f(t),  surface  tractions  t(t)  on  the  boundary  Fj  and  prescribed  dis¬ 
placement  fields  on  the  boundary  Fa.  In  real  heterogeneous  materials,  dimensions  of  the 
RME  of  characteristic  length  I  are  typically  very  small  in  comparison  with  the  dimensions  of 
the  body  of  characteristic  length  L.  The  ratio  of  these  microscopic  and  macroscopic  scales  //L 
is  represented  by  a  very  small  positive  number  e.  When  subjected  to  structural  level  loads 
and  displacements,  the  resulting  evolutionary  variables  e.g.  deformation  and  stresses,  vary 
from  point  to  point  at  the  macroscopic  scale  x.  Furthermore,  a  high  level  of  heterogeneity 
in  the  microstructure  causes  a  rapid  variation  of  these  variables  in  a  small  neighborhood  e 
of  the  macroscopic  point  x.  This  corresponds  to  a  microscopic  scale  -  and  consequently, 
all  variables  are  assumed  to  exhibit  dependence  on  both  length  scales  i.e.  ^ (x,  — ). 

The  superscript  e  denotes  association  of  the  function  with  the  two  length  scales.  In  this 
notation,  f)"  corresponds  to  a  connected  domain  that  extends  the  structural  domain  to  its 
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(c) 

Figure  4.1:  A  heterogeneous  structure  with  various  levels  (a)  The  global  structure  and  compu¬ 
tational  geometry  (b)  Representative  material  elements  (RME)  at  a  point  and  the  corresponding 
VCFE  model  and,  (c)  A  basic  structural  element  (BSE)  represented  by  a  Voronoi  Cell. 


50 


microstructure.  Mathematically  speaking. 


fi'  =  {x  €  Q  :  0(^)  =  1}  (4.1) 

in  which  0(y)  =  1  when  y  lies  in  the  microscopic  RME.  In  most  of  the  work  on  homogeniza¬ 
tion  theory  [22,  18,  20,  25],  a  periodic  repetition  of  the  microstructure  about  a  macroscopic 
point  X  has  been  assumed,  thereby  making  the  dependence  of  the  function  on  y  =  (-),  peri¬ 
odic.  This  characteristic  is  often  termed  as  Y-periodicity,  where  Y  corresponds  to  a  RME. 
The  instantaneous  tangent  modulus  tensor  relating  the  stress  rate  (increment)  to  the 

strain  rate  (increment),  and  the  corresponding  instantaneous  compliance  tensor  Sljki  in  the 
connected  domain  are  expressed  as: 

=  Ep«(x,y)  in  (4.2) 

^  E-]liix,y)  in  (4.3) 

It  is  assumed  that  the  stress,  strain  and  displacement  fields  satisfy  the  rate  or  incremental 
forms  of  the  equilibrium  equation,  kinematic  relation,  and  constitutive  relations  given  as; 


2  V(9a;^  dxl) 


in  9^ 


(4.5) 


=  Etjkkli  in  99  (4.6) 

where  u\  =  ■U;(x,  y)  is  a  Y-periodic  rate  of  displacement  field  in  y.  Furthermore  the  boundary 
conditions  are  assumed  to  satisfy  the  following  equations  on  the  prescribed  traction  and 
displacement  boundaries  respectively. 

(T-yUj  =  ii  on  Tt  (4.7) 

u\  =  u,  on  r„  (4.8) 

where  n  is  the  unit  normal  to  the  boundary.  Since  small  deformation  is  assumed  in  this 
analysis,  the  normal  vector  does  not  change  significantly  with  load.  In  homogenization 
theory,  the  Y-periodic  displacement  rate  or  increment  field  is  approximated  by  an  asymptotic 
expansion  with  respect  to  parameter  e: 


u^(x)  =  u°(x,y)  +  eui(x,y)  +  eM^(x,y)  +  ---,  y=^  (4.9) 


Noting  that  the  spatial  x^  derivative  of  any  function  depends  on  the  two  length  scales  and 


IS  given  as: 


(4.10) 
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the  strain  rate  tensor  may  be  expressed  as: 


1 du\ 
t  dyi  dyi 


du\  dill 
dxi  dyi 


(4.11) 


Substituting  this  in  the  constitutive  relation  (4.6),  the  stress  rate  ajj  can  be  expanded  as 


where 


-dj  +  dj.  +  eal  +  e^afj 


•0  _  pt 

ijkl 

(J.  .  _  i"  o 


(4.12) 


^^‘\dxi  dyi 
,  /dill  dill 
[dxi  ^  dyi 


(4.13) 


Putting  the  expansion  of  (4.12)  in  the  rate  form  of  equilibrium  equation  (4.4),  and  setting 
each  coefficient  of  e*,  {i  =  —1,0, 1,2, ...)  to  zero,  results  in  the  following  set  of  equations. 


dafj 

=  0 


dajj  da^^  _ 
dyj  dxj 


dyj  dxj 


+  /«  —  0 


(4.14) 


If  X  and  y  are  considered  as  independent  variables,  equations  (4.14)  form  a  recursive  sys¬ 
tem  of  differential  equations  for  the  functions  uj,uf...  parameterized  by  x.  An  important 
result  that  is  used  to  solve  the  system  of  equations  (4.14)  is  given  as  follows  [24].  The  equation 


(4.15) 


for  a  Y-  periodic  function  $  =  $(x,  y)  has  a  unique  solution  if  the  mean  value  of  F,  defined 

by 

<F>=j^/^Fdy  =  0  (4.16) 

where  |Y|  denotes  the  volume  of  the  RME.  Equations  (4.14)  and  (4.13)  together  with  results 
(4.16)  leads  to  the  trivial  value  for  afj,  and  therefore  establishes  that  u°  is  only  a  function 
of  X  as  shown  in  [25],  i.e. 

dfj  =  0,  and  =  ^?(x)  (4.17) 
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4.2.1  Microscopic  Equilibrium  Equation  and  Homogenized  Constitutive  Rela¬ 
tion 


Substituting  (4.17)  into  the  second  of  equation  (4.14)  leads  to  the  Y-  domain  equilibrium 
equation 


dyj 


0 


(4.18) 


From  equations  (4.6),  (4.11)  and  (4.12),  by  neglecting  the  terms  associated  with  e  or  higher, 
the  constitutive  relation  in  Y  is  expressed  as 


(4.19) 


where 


'kl 


+  ^kl  — 


du\  dill 

dxt  dyi 


(4.20) 


Here  eh  is  the  local  or  microstructural  strain  rate  tensor,  for  which  Cki  =  is  an  averaged 

du^ 

macroscopic  part,  and  eh  —  is  denoted  as  a  fluctuating  strain  rate  tensor  (Suquet  [32]). 
With  the  equilibriated  microscopic  stress  rate  and  displacement  rate  ul  fields  completely 
determined  by  equations  (4.18),  (4.19)  and  (4.20),  a  local  structure  tensor  denoted  by 
is  calculated  by  applying  a  unit  components  of  a  macroscopic  strain  rate  tensor.  The  local 
structure  tensor  is  used  to  establish  a  relation  between  the  local  microscopic  strain  rate 
eh,  and  the  macroscopic  strain  rate  e*;  as: 


^kl  ^pm^pm 


(4.21) 


Due  to  linearity  of  the  rate(incremental)  problem,  a\-  and  ii]  can  be  expressed  in  the  forms 


where 


dVi 


dxi' 


u] 


x{y)i 


kl 


dxi 


{microscopic  equilibrium) 


(4.22) 

(4.23) 


In  equation  (4.22),  is  a  Y-antiperiodic  function  and  is  a  Y-periodic  function  represent¬ 
ing  characteristic  modes  of  the  RME.  Substituting  equations  (4.22)  and  (4.21)  in  equation 
(4.13)  yields  the  microscopic  constitutive  relations  as: 


<^q(y)  =  (4.24) 

The  local  structure  tensor  is  related  to  the  Y-periodic  function  Xi‘  as: 

Krn  =  K  +  ^1  (4.25) 
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where  T-f  is  a  fourth  order  identity  tensor  expressed  as: 


Tt'  =  +  6ui,i,) 


(4.26) 


The  set  of  equations  (4.22)  and  (4.25)  determine  the  vector  x(y)fc  to  within  an  additive  con¬ 
stant.  The  mean  of  equation  (4.25)  yields  the  homogenized  elastic-plastic  tangent  modulus, 
for  use  in  the  macroscopic  analysis,  in  the  form 


^ijkl  ^  ^ij  ^ 


= 


r 

Iy 


(4.27) 


4.2.2  Macroscopic  Equations 

Taking  the  mean  of  the  third  equation  in  (4.14)  on  Y  and  applying  the  condition  (4.16)  to 
the  first  term  leads  to  an  averaged  form  of  the  global  equilibrium  equation  as: 


d  <  >  ;  - 

+/.  =  0 

OXj 


ITl  ^structure 


(4.28) 


Note  that  the  above  equation  (4.28)  is  now  valid  for  the  macroscopic  domain  ^structure-  Thus 
in  the  macroscopic  domain,  the  averaged  stress  rate  S  =<  >  and  displacement  rate  u° 

fields  are  the  solutions  to  the  elastic-plastic  problem  delineated  by  the  equations: 

—  f  ■  O 

—  Ji  ^  ^structure 


Stj  —  E-jj^iCki 

t,ijnj  =  ii  on  Ft 

•  0  —  "p 

u,  =  Ui  on  I  u 


(4.29) 


An  incremental  small  deformation  analysis  is  pursued  with  the  equations  (4.29)  at  the 
macroscopic  level,  and  by  using  the  Voronoi  cell  finite  element  model  (VCFEM)  for  solving 
the  microscopic  problem.  Developments  of  the  incremental  VCFEM  model  for  heterogeneous 
microstructures  are  presented  in  the  next  section. 


4.3  Microstructural  Analysis  with  VCFEM 

Voronoi  cells,  resulting  from  Dirichlet  tessellation  of  a  heterogeneous  microstructure,  make 
rather  unconventionaJ  elements,  due  to  the  arbitrariness  in  the  number  of  edges.  The  Voronoi 
Cell  finite  element  model  developed  by  Ghosh  and  coworkers  [26,  27,  28,  29]  avoids  difficulties 
of  conventional  displacement  based  FEM  formulations  by  invoking  the  assumed  stress  hybrid 
method  ,introduced  by  Pian  [33].  In  this  formulation,  independent  assumptions  are  made 
on  an  equilibriated  stress  field  in  the  interior  of  each  element  and  a  compatible  displacement 
field  on  the  element  boundary.  Small  deformation  elastic-plastic  analysis  of  materials  with 
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embedded  second  phase  has  shown  significant  promise  with  respect  to  efficiency  and  accuracy 
[26,  27].  Details  of  this  development  based  on  the  hybrid  formulation,  originally  proposed 
by  Atluri  and  coworkers  [34,  35]),  are  presented  in  [26,  27].  In  this  section,  a  brief  account 
of  this  formulation  is  presented  for  completeness. 

4.3.1  Variational  Principles  in  VCFEM 

Consider  a  typical  representative  material  element  (RME)  Y  tessellated  into  N  Voronoi  cells, 
as  shown  in  figure  4.1(b).  This  is  based  on  the  location,  shape  and  size  of  N  heterogeneities 
as  explained  in  [30].  The  matrix  phase  in  each  Voronoi  cell  Yg  is  denoted  by  Ym  and  the 
heterogeneity  (void  or  inclusion)is  denoted  by  K.  The  matrix-heterogeneity  interface  dYc 
has  an  outward  normal  n'",  while  n^  is  the  outward  normal  to  the  element  boundary  dYg. 

In  this  presentation,  n  incremental  finite  element  formulation  is  invoked  to  account  for  rate 
independent  plasticity.  At  the  beginning  of  the  p-th  increment,  let  cr  be  an  equilibriated 
stress  field  with  a  strain  field  e[<r^load  history),  and  u  be  a  compatible  displacement  field 
on  the  element  boundary.  Also  let  A<t  correspond  to  an  equilibriated  stress  increment  in 
Ye,  Au  to  a  compatible  displacement  increment  on  dYg,  and  At  to  a  traction  increment  on 
the  traction  boundary  Ftm-  The  incremental  problem  is  solved  by  using  a  two  field  assumed 
stress  hybrid  variational  principle,  derived  from  an  element  energy  functional  as: 

ne(A<T,  Au)  -  -  /  AB{o-,  Aa)  dY  -  [  e  :  Aa  dY 

JYe  JYe 

+  [  {<T  +  A(r)  ■  ■  {u  +  Au)  dY  -  f  (t  +  At)  •  (u -|- Au)  dF 

•tdYe  JTtm 

-  (<7”^  +  A^t”*  -  -  A<t")  •  n"  •  (u'  +  Au')  dY  (4.30) 

where  Au'  is  the  displacement  of  the  interface  and  AB  is  the  increment  in  element  compli¬ 
mentary  energy.  Superscripts  m  and  c  represent  respectively  the  matrix  and  second  phase 
parts  of  the  Voronoi  cell  element.  The  energy  functional  for  the  entire  domain  is  obtained 
by  adding  each  element  functional  as 

N 

n  =  (4.31) 

e=l 

The  first  variation  of  lie  with  respect  to  the  stress  increments  A<t,  results  in  the  kinematic 
relations  as  the  Euler  equation, 

VAu  =  Ae  in  Ye  (4.32) 

while  the  first  variation  of  FI  with  respect  to  boundary  displacement  increments  Au  yields 
traction  conditions  as  Euler  equations, 

{cr  +  Act)  ■  n®"''  =  — (o-  -|-  A<t)  ■  on  cell  boundary  F^  (Interelement  traction  reciprocity) 
(<7  -|-  A<t)  -11®=  t  -b  At  on  traction  boundary  F^^ 

+  A<r‘')  •  n‘'  =  (cr'^  +  A<r™)  •  on  interface  dYc  (Interface  traction  reciprocity)  (4.33) 


Equilibriated  stress  increments  A<t,  constitutive  relations,  along  with  the  incremented  form 
of  the  energy  functional  completely  define  the  microstructural  problem  in  the  p-th  increment. 


Element  formulations  and  assumptions 

In  the  Voronoi  cell  finite  element  model  (VCFEM)  formulation,  independent  assumptions 
on  stress  increments  A<r  are  made  in  the  matrix  and  heterogeneity  phases  to  accommodate 
stress  jumps  across  the  interface.  Use  of  stress  functions  $(a;,y)  is  a  convenient  way  of 
deriving  stress  increments  in  two-dimensional  analysis.  Different  expressions  may  be  assumed 
for  $  in  the  matrix  and  inclusion  phases.  In  general,  these  can  be  arbitrary  functions  of 
location,  yielding  stress  increments  in  the  form 

{  A<r-  }=  [p-(a:,t/)]{  A/3”^  } 

{  A<t^  }  =  [  P‘=(x,j/)  ]  {  A/3^  }  (4.34) 

where  { A/3}’s  correspond  to  a  set  of  yet  undetermined  stress  coefficients  and  [P]  is  a  matrix 
of  interpolation  functions.  Compatible  displacement  increments  on  the  element  boundary 
dYe  as  well  as  on  the  interface  51^,  are  generated  by  interpolation  in  terms  of  generalized 
nodal  values.  The  displacement  increments  on  the  element  boundary  and  interface  may  then 
be  written  as, 

{Au}  =  [Ll{Aq} 

{Au'}  =  [L"]{Aq'}  (4.35) 

where  {Aq}  and  {Aq'}  are  generalized  displacement  increment  vectors  and  [L]  is  an  inter¬ 
polation  matrix.  Substituting  element  approximations  for  stresses  (4.34)  and  displacements 
(4.35),  in  the  energy  functional  (4.30),  and  setting  the  first  variations  with  respect  to  the 
stress  coefficients  AyS”*  and  A/3‘'  respectively  to  zero,  results  in  the  following  two  weak  forms 
of  the  kinematic  relations  (4.32), 

/  [P"]^{e -h  Ae}dF  =  /  [P’"]^H[L^]dy  {Aq}  -  /  [P”^f  [n‘^][L‘=]dF  {Aq'} 

JYm  *JdYc,  JdYc 

[  [P^]^{e  + Ae}dF  =  /  [P"]^[n‘=][L"]dy  {Aq'}  (4.36) 

JYc  JdYc 

Setting  the  first  variation  of  the  total  energy  functional  (4.31)  with  respect  to  Aq  and 
Aq'  to  zero,  results  in  the  weak  form  of  the  traction  reciprocity  conditions  as 

r  [P”lrfy  0  1  /  +  A,3“  \  _ 

s[-/3nW[n'r[P"l<ir  V.|V)qn'|^[P1rfy  J  \  +  f- 


/r,„[L'F{t  +  AtjdV 

m 


(4.37) 


For  an  elastic-plastic  material,  the  strain  increments  Ae  in  equation  (4.36)  are  non-linear 
functions  of  the  current  state  of  stress  cr  as  well  as  of  their  increments  A<r.  The  non-linear 
finite  element  equations  (4.36)  and  (4.37)  are  solved  for  the  stress  parameters  (A/B”*,  A/S”^) 
and  the  nodal  displacement  increments  (Aq,  Aq')  at  the  p-th  increment. 


4.4  Elastic-plastic  Homogenization  with  VCFEM 


In  this  section,  the  Voronoi  Cell  finite  element  model  is  formulated  to  be  applied  in  con¬ 
junction  with  the  homogenization  method  for  coupling  global  and  local  analyses.  VCFEM 
is  used  to  model  an  arbitrary  microstructural  RME,  and  consequently,  Y  represents  a  RME 
with  a  boundary  5Y.  In  an  incremental  formulation,  the  equilibriated  microscopic  stress 
increment  corresponds  to  A<t^(=  A<t^)  in  equation  (4.22)  and  the  microstructural  strain 
increments  are  designated  as  Ae^  in  equation  (4.20).  Similarly,  the  increments  in  micro¬ 
scopic  displacements  on  the  cell  boundaries  dYg  are  identified  with  Aw^  in  equation  (4.22) 
and  those  on  the  interface  are  denoted  by  Att^'.  In  the  absence  of  traction  boundaries  due 
to  periodicity  conditions,  the  incremental  energy  functional  for  each  Voronoi  cell  element  in 
equation  (4.30)  is  modified  for  the  homogenization  process  as: 


~Ir.  („‘  +  A^)  n‘  ddY 

L  (“t  +  +  f  (««  +  Ae.v)Acr'p(K.38) 

JYe 


where  5/^^  is  an  instantaneous  elatic-plastic  compliance  tensor.  The  last  term  in  equation 
(4.38)  incorporates  the  effect  of  macroscopic  strains  in  the  microstructure.  The  station¬ 
ary  condition  of  He  with  respect  to  stress  increment  Acr^^  yields,  as  Euler’s  equations,  the 
incremental  form  of  kinematic  relations  (4.20) 


+  Ae,j  —  e,j 


Ae,-,'  + 


d(ul  +  An-) 


Furthermore,  stationarity  of  the  total  energy  functional  11  =  Re  with  respect  to 
displacement  increments  Auj  and  uj'  respectively,  result  in  the  inter-element  and  interface 
traction  reciprocity  conditions: 


(af^  +  A(7^)-<  =  -K  +  A<T'^.)-^r  (4.40) 

(cr';  +  Act';)  •  n]  =  (cr'™  -f  Act'™)  •  n'  on  dY  (4.41) 

where  superscript  +  and  -  denote  values  on  opposite  sides  of  the  inter  element  boundary 
dYg.  The  three  Euler’s  equations  together  with  (a)  the  assumed  equilibriated  stress  fields 
satisfying  Acr^  -  =  0  in  (b)  assumed  compatible  displacement  increment  fields  in  dY^  and 
dYc,  and  the  instantaneous  constitutive  relation  Ae\-  =  S''j^(A(t^;  describe  the  incremental 
boundary  value  problem  for  the  microstructure. 
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The  microscopic  VCFEM  module  is  executed  for  two  purposes  in  each  increment  of  the 
macroscopic  module.  The  first  is  to  evaluate  the  microscopic  stress  increments  Act'  from 
given  values  of  the  macroscopic  strain  e  at  the  beginning  of  the  step,  and  its  increment  Ae. 
The  second  is  to  calculate  the  instantaneous  homogenized  tangent  modulus  at  the  end 
of  the  increment  in  the  macroscopic  module. 

4.4.1  Calculation  of  Microscopic  Stresses 

This  step  involves  the  iterative  solution  of  modified  forms  of  the  kinematic  relations  (4.36) 
and  the  traction  reciprocity  conditions  (4.37)  to  yield  the  incremental  stress  parameters 
A/3  and  the  nodal  displacement  increments  Aq  and  Aq^  The  weak  form  of  the  element 
kinematic  relation  (4.39)  is  obtained  by  modifying  (4.36)  for  homogenization  in  the  multiple 
scale  model  as: 


/  [P™F{e^  +  Ae^  -  e  -  Ae}d  y  =  /  [n^][L^]dr  {q  +  Aq}e 

-  Yc  JdYe 

-  [  [PY[n1[Lldy  {q'  + Aq'},  (4.42) 

JdYc 

I  [Pl^{e'  +  Ae'  -  e  -  Ae}d  Y  =  f  [P']^[n‘=][L']dy  {q'  +  Aq'},  (4.43) 

As  mentioned  in  section  3.2,  an  iterative  solution  process  is  implemented  for  evaluating 
stresses  from  given  values  of  the  nodal  displacements.  In  the  i-th  iteration,  the  total  strain 
increment  {Ae'}  is  linearized  with  respect  to  stress  increment  {Act'}  in  the  form: 

{Ae'}  =  {Ae'(o-',  Ao-')}’  +  [S']  :  {do-'}‘  (4.44) 

where  the  correction  {d<r'}‘  to  stress  increment  is  expressed  as, 

{d<r'}‘  =  [P]{d/3y  in  F,  (4.45) 

Substitution  in  (4.42)  and  (4.43)  leads  to  the  linearized  kinematic  equations  (3.22)  in  con¬ 
junction  with  the  homogenization  procedure  as: 

■  H,,  0  1  /  d/3’"  y  _  r  G,  1  r  q  Aq  1  _  r  /K._KjP’"]^{e'  +  Ae'’  -  e  -  Ae}d  F  1 

0  H,  J\d/3^j"[0  G„  J  \  q'-H  Aq'  j  \  /i.jP']^{e' +  Ae'*  -  e  -  Ae}d  F  j 

(4.46) 

where  components  of  the  matrices  are  explained  in  equation  (3.22).  The  traction  reciprocity 
equations  (3.23)  are  iteratively  solved  for  nodal  displacements  in  the  same  way  a.s  mentioned 
in  section  3.3.  In  the  j-th  iteration,  this  corresponds  to  the  solution  of  the  linearized  global 
traction  reciprocity  equation: 

£[GnHl-[G|  I  I’  = 

_  f  [  /anlI'1VFlP“F<iV  0  1  f /J”  +  1 ' 
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4.4.2  Calculation  of  Homogenized  Tangent  Modulus 

Evaluation  of  the  homogenized  tangent  modulus  tensor  is  performed  at  the  end  of  an 
increment  in  the  global  analysis  (ABAQUS)  after  converged  values  of  microscopic  variables 
have  been  obtained.  From  the  basic  definition  in  equation  (4.27),  components  of  are  ob¬ 
tained  by  averaging  the  true  stress  increments  in  response  to  unit  increments  in  components 
of  the  strain  tensor  e.  In  two  dimensional  analysis  the  individual  components  correspond  to 
611,622  and  612,  or  in  contracted  notations  61,62  and  §3  i.e.  e,-,  i  =  1,2,3.  This  is  obtained 
first  by  consistent  linearization  of  the  element  kinematic  relation  in  equation  (3.22)  with 
respect  to  each  component  of  the  macroscopic  strain  increment  6;  to  yield: 


0  1  /  Sr 
0  H,  J  \  8^^ 


Ge  —Gem  I  Sq 
0  Gee  \  \  Sq' 


Sr.-Y,{'P''V  {Si.}dY 

/yJP=]^{te,}dF 


4.48) 


In  equation  (4.48),  it  is  assumed  that  the  stress  coefficients  /3  and  generalized  nodal  dis¬ 
placements  q  are  perturbed  about  their  respective  equilibriated  and  compatible  values,  cor¬ 
responding  to  a  prescribed  perturbation  in  e,-.  The  nodal  displacements  in  equation  (4.48) 
are  solved  by  linearization  of  the  global  traction  reciprocity  condition  (4.37)  with  respect  to 
Oj-,  given  as: 


E(gF  ((h)-‘[g|  1 I 


[«]-■ 


{8-e,}dY 


=  0  (4.49) 


The  solved  values  of  {dq},-  yield  the  characteristic  modes  of  deformation  as  plotted  in  figures 
(4.14),  (4.15)  and  (4.16).  These  are  subsequently  substituted  in  the  local  compatibility 
equation  (4.48)  to  solve  for  element  stress  coefficients  {^/3},.  It  should  be  noted  the  tangent 


modulus  is  evaluated  from  unit  values  of  {(^ej,  i.e.  {(5ei}  = 

f  0  1 


1 1 

0  1 

0 

^  1  {^62}  —  •' 

1 

0  J 

[o  1 

,{^63}  = 


0 

1 


Components  of  the  elatic-plastic  homogenized  tangent  modulus  tensor  Ef^j^  are  then 


obtained  from  equation  (4.27),  by  averaging  [P]{d/3},  for  each  component  of  the  macroscopic 
strain,  as 


[Eq  =  i 


E/jPlWh  E/lPlP/sh  i:/[P|W3 

e=l“'E<=  e=mP'=  e=l 


(4.50) 


4.5  Numerical  Implementation 

A  number  of  numerical  details  have  to  be  considered  in  the  construction  of  a  multiple  scale 
computational  model  using  the  microstructural  VCFEM  model.  In  this  section,  a  few  salient 
features  are  discussed. 
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4.5.1  Incorporation  in  the  Macroscopic  Analysis  Module 

The  Voronoi  cell  finite  element  module  is  incorporated  in  a  macroscopic  analysis  module 
with  the  interface  being  created  by  the  homogenization  procedure.  In  this  work,  the  general 
purpose  commercial  code  ABAQUS  has  been  chosen  to  serve  as  the  macroscopic  analy¬ 
sis  program.  Macroscopic  ABAQUS  models  are  developed  with  two  dimensional  elements. 
The  material  constitutive  relation  at  each  integration  point  of  ABAQUS  elements  is  input 
through  the  homogenization  process  by  using  results  from  the  microscopic  VCFEM  analy¬ 
sis.  This  interface  between  ABAQUS  and  VCFEM  is  created  through  UMAT,  intended  to 
incorporate  user  specified  constitutive  models  in  ABAQUS.  The  analysis  code  resulting  from 
this  macro-micro  coupling  is  termed  as  VCFEM-HOMO.  QUAD4  elements  with  one-point 
reduced  integration  and  hourglass  control  are  selected  in  ABAQUS  analysis.  ABAQUS  and 
VCFEM  codes  used  in  this  analysis  implement  implicit  time  integration  schemes,  and  hence 
require  iterations  at  both  levels  for  convergence. 

Within  each  iteration  loop  of  macroscopic  analysis  in  the  p-th  increment,  microscopic 
state  variables  are  computed  in  VCFEM  using  given  values  of  the  macroscopic  strains  at  the 
beginning  of  the  increment  and  its  increment,  as  shown  in  the  flow  chart  of  flgure  (4.2).  The 
microscopic  stresses  are  averaged  to  yield  macroscopic  stresses.  Upon  convergence  in  global 
and  local  iterations,  the  microscopic  VCFEM  is  invoked  again  to  evaluate  the  homogenized 
elastic-plastic  tangent  modulus  by  applying  unit  components  of  macroscopic  strain.  All 
the  history  dependent  microstructural  variables,  e.g  stresses,  strains  and  plastic  strains  are 
retained  and  updated  during  the  incremental  process. 

4.5.2  Implementing  Periodicity  Boundary  Condition 

An  essential  step  in  computing  the  homogenized  material  properties  is  to  ensure  repeatability 
boundary  conditions  on  the  RME.  For  a  square  or  rectangular  RME,  identical  displacement 
functions  must  be  specified  for  corresponding  nodes  (equidistant  from  a  coordinate  axis) 
on  opposite  edges.  Implementing  the  repeatability  boundary  condition  for  a  regular  finite 
element  mesh  is  straightforward,  since  a  uniform  mesh  can  be  generated  to  have  the  cor¬ 
respondence  between  boundary  nodes  on  opposite  faces.  However  boundary  nodes  of  the 
Voronoi  mesh  generated  by  Dirichlet  tessellation,  are  in  general  quite  arbitrary  and  such 
a  correspondence  cannot  be  easily  established.  A  method  which  involves  representation  of 
nodal  boundary  displacements  by  a  suitable  polynomial  function  is  implemented  for  enacting 
the  repeatability  conditions.  In  this  method,  a  (^  —  l)th  order  polynomial  is  chosen  for  the 
displacements,  where  q  corresponds  to  the  highest  number  of  boundary  nodes  between  the 
two  opposite  faces.  That  is,  if  one  face  has  5  nodes  while  the  opposite  side  consists  of  6 
nodes,  a  5-th  order  polynomial  function  is  chosen.  The  edge  nodal  displacements  are  then 
written  as: 


Ui  —  Co  4“  -fi  02^1  T  (^33^1  4"  ci4X^  4" 

=  60  4-  biVi  4-  62J/1  4-  &32/i  4-  64J/1  +  b^yl 
U2  =  cto  4-  0,1X2  +  02x1  -I-  03X2  4-  04X2  4-  05X\ 


60 


Figure  4.2:  Flow  chart  of  microscopic  analysis 
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V2  =  bo  +  612/2  +  ^22/2  +  ^32/2  +  ^42/2  +  ^52/2 


(4.51) 


where  ui,  Vi,U2,  V2  •  •  •  are  nodal  displacements  and  xi,  i/i,  X2, 2/2  ••  •  are  boundary  coordinates. 
These  lead  to  displacement  constraints  in  the  matrix  equations  prior  to  solving. 


4.6  Numerical  Examples 

Numerical  examples  conducted  with  the  multiple  scale  homogenization  module  VCFEM- 
HOMO,  are  divided  into  two  categories.  The  first  set  of  examples  are  intended  to  validate 
the  effectiveness  of  the  asymptotic  homogenization  in  conjunction  with  Voronoi  cell  FEM 
formulation  for  heterogeneous  materials  with  both  inclusions  and  voids.  This  is  accomplished 
through  comparison  of  VCFEM-HOMO  results  in  simple  tests,  with  (a)  those  generated  by 
Unit  cell  models  using  conventional  finite  element  codes  such  as  ABAQUS,  (b)  experimen¬ 
tal/analytical  results,  and  (c)  effective  continuum  models.  The  comparisons  also  help  in 
identifying  some  of  the  shortcomings  of  alternative  approaches,  and  where  coupled  multiple 
scale  analyses  is  desirable.  In  the  second  set  of  examples,  the  VCFEM-HOMO  code  is  used 
to  solve  more  complex  multiple  scale  problems  with  various  microstructural  morphologies. 
Effects  of  the  microsctructure  in  the  evolution  of  macroscopic  and  microscopic  variables  are 
investigated.  The  material  in  all  the  ensuing  examples  is  assumed  to  be  a  boron-aluminum 
composite  or  aluminum  with  voids,  unless  otherwise  mentioned.  The  boron  inclusion  is 
assumed  to  be  elastic  while  the  aluminum  matrix  is  an  elastic-plastic  material  with  the  fol¬ 
lowing  properties. 

Material  Properties: 

Boron  fiber 

Young’s  Modulus  {Ec):  344.5GPa  Poisson  Ratio  (uc):  0.26 

Aluminum  matrix 

Young’s  Modulus  {Em)'-  68.9  GPa  Poisson  Ratio  {I'm)-  0.32 

Initial  Yield  Stress  (5^):  94  MPa  Post  Yield  flow  rule  ;ee,„  = 

The  microstructural  VCFEM  analysis  uses  a  48  term  stress  function  in  equation(3.29) 
(p  -f  ^  =  2. .4  and  i  =  1..3)  for  the  matrix  of  porous  materials,  and  a  34  term  stress  function 
(p  -|-  ^  =  2  and  i  =  1..3)  in  the  matrix  of  composite  materials.  The  stress  field  in  the  in¬ 
clusion  for  composite  materials  is  generated  with  a  25  term  polynomial  stress  function  (i.e. 
^poiy'  P  +  9  =  2.. 6)). It  should  be  noted  that  elastic  problems  with  homogenization  and 
VCFEM  have  been  successfully  solved  in  [28]. 

4.6.1  Validation  of  the  Homogenization  Model 
Numerical  Unit  cell  models 

In  this  example,  results  of  VCFEM-HOMO  are  compared  with  predictions  of  numerical  unit 
cell  models  for  plane  strain  elastic-plastic  analysis.  Two  microstructural  arrangements  viz. 
square  edge  distribution  (SED)  and  hexagonal  distribution  (HD)  as  shown  in  figure  4.3,  are 
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considered.  For  VCFEM-HOMO  analysis,  the  macroscopic  ABAQUS  mesh  consists  of  a  sin¬ 
gle  QUAD4  element,  and  a  single  Voronoi  cell  element  is  required  for  microstructural  analysis 
as  shown  in  figure  4.3(a)  and  (c).  Unit  cell  analyses  are  conducted  at  the  microstructural 
level  only,  with  the  representative  material  element  (RME)  modeled  by  QUAD4  elements 
in  ABAQUS,  illustrated  in  figure  4.3(b)  and  (d).  The  inclusion  volume  fraction  is  55%. 
The  periodic  construct  of  unit  cell  is  achieved  through  restrictive  repeatability  conditions 
enforced  on  edges  x  =  L  and  y  =  L.  This  constraint  condition  requires  rectangular  cells 
to  deform  into  rectangular  shapes,  i.e.  straight  line  edges  {x  =  L,y  =  L)  move  uniformly 
as  straight  lines.  For  the  SED  and  HD  arrangements,  the  ABAQUS  mesh  consists  of  1440 
QUAD4  and  1482  elements  respectively. 

Comparisons  are  made  for  macroscopic  as  well  as  microscopic  response  by  the  two  ap¬ 
proaches.  For  the  unit  cell  model,  macroscopic  variables  denoted  with  an  overbar,  are  ob¬ 
tained  by  volume  averaging  the  microscopic  response.  For  example,  macroscopic  stresses  & 
are  evaluated  as: 


Jo  crdQ 

^  =  27  =  - 

fo  dn 

Figures  4.4  and  4.5  show  the  macroscopic  stress-strain  relations  and  the  microscopic  stress 
distribution  for  the  square  edge  and  hexagonal  packings  respectively  when  subjected  to 
uniform  stretching.  Figures  4.4(a)  and  4.5(a)  show  excellent  agreement  between  VCFEM- 
HOMO  and  Unit  Cell  model  in  the  macroscopic  responses  at  all  load  stages.  Figures  4.4(b) 
and  4.5(b)  show  the  microscopic  stress  distribution  in  the  direction  of  applied  strain  along 
sections  at  y  =  ^  and  t/  =  0  in  figures  4.3(b)  and  (d)  respectively.  Though  the  VCFEM- 
HOMO  results  for  both  packings  compare  well  with  unit  cell  predictions,  better  agreement 
is  seen  with  HD  packing  due  to  the  refined  unit  cell  mesh  along  the  section  considered. 

VCFEM-HOMO  code  is  next  utilized  to  investigate  the  effect  of  microstructural  mor¬ 
phology  on  the  material  response.  Four  different  arrangements  are  considered  with  circular 
inclusions  of  55%  volume  fraction,  viz.  (a)  square  edge  distribution  (figure  4.3(a)),  (b)  hexag¬ 
onal  distribution  (figure  4.3(c)),  (c)  square  diagonal  distribution  (45°  rotation  of  figure  4.3(a)) 
and  (d)  random  distribution  with  29  inclusions  (figure  4.6)  The  ABAQUS  macroscopic  mesh 
consists  of  one  plane  strain  element,  subjected  to  simple  tension  and  simple  shear  boundary 
conditions  to  maximum  effective  strains  of  0.008.  The  effective  macroscopic  stress-strain 
behavior  is  depicted  in  figure  4.7.  On  account  of  the  rotational  equivalence  of  RMEs,  the 
behavior  of  the  square  edge  and  square  diagonal  distributions  are  exactly  reversed  for  the  two 
loading  conditions.  The  results  indicate  the  strong  anisotropy  in  overall  behavior  of  these 
two  RME’s.  The  random  and  hexagonal  packings  on  the  other  hand  exhibit  similar  response 
in  tension  and  shear.  This  establishes  the  low  directional  dependency  of  these  RME’s,  im¬ 
plying  near  isotropic  behavior.  The  stress-strain  curves  for  the  latter  two  distributions  lie 
within  the  envelopes  of  the  SE  and  SD  distributions,  and  hence  exhibit  more  overall  ductility. 
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Figure  4.3:  Computational  model  of  a  composite  material  with  a  circular  inclusion  of  volume 
fraction  55%  (a)  Macroscopic  and  microscopic  models  for  Square  edge  distribution  by  VCFEM- 
HOMO,  (b)  Microstructural  ABAQUS  Unit  Cell  model  for  Square  edge  distribution  (c)  Macroscopic 
and  microscopic  models  for  Hexagonal  distribution  by  VCFEM-HOMO,  and  (d)  Microstructural 
ABAQUS  Unit  Cell  model  for  Hexagonal  distribution. 


strain 


(b) 


Figure  4.4:  Comparison  of  macroscopic  and  microscopic  response  for  Square  edge  packing  composite 
between  VCFEM-HOMO  and  Unit  Cell  models  (a)  macroscopic  stress-strain  response  along  the 
load  direction  and  (b)  microstructural  stress  distribution  along  a  section  through  the  inclusion  at 

X  =  0.0 
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(a) 


y/L1 

(b) 

Figure  4.5:  Comparison  of  macroscopic  and  microscopic  response  for  Hexagon  packing  composite 
between  VCFEM-HOMO  and  Unit  Cell  models  (a)  macroscopic  stress-strain  response  along  the 
load  direction  and  (b)  microstructural  stress  distribution  along  a  section  through  the  inclusion 
a:  =  0 
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Effective  continuum  models 

This  set  of  examples  examines  some  of  the  differences  between  effective  continuum  models 
and  coupled  macro-micro  analysis  for  simple  loadings.  Furthermore,  the  VCFEM-HOMO 
code  is  also  tested  against  experimental  and  analytical  results. 

In  the  first  example,  overall  stress-strain  plots  are  first  generated  by  uniaxial  stretching 
of  VCFEM  models  in  figures  4.3(a)  and  4.6  under  plane  stress  conditions.  This  is  then 
input  as  the  homogenized  material  law  in  ABAQUS  with  isotropic  hardening  and  associ¬ 
ated  flow  rule.  The  ABAQUS  mesh  with  this  effective  continuum  model  consists  of  one 
QUAD4  element.  Likewise,  the  macroscopic  ABAQUS  mesh  for  VCFEM-HOMO  consists 
of  one  QUAD4  element  and  the  microscopic  models  are  illustrated  in  figures  4.3(a)  and  4.6. 
Macroscopic  longitudinal  stress-strain  relations,  and  transverse-longitudinal  strain  relations 
by  the  two  approaches  are  plotted  in  figures  4.8  and  4.9.  Comparisons  are  made  for  three 
different  microstructures,  viz.  (a)  square  edge  packing  with  10%  inclusion  volume  fraction, 
(b)  square  edge  packing  with  55%  inclusion  volume  fraction,  and  (c)  random  packing  with 
55%  inclusion  volume  fraction.  The  stress  and  strain  plot  in  the  loading  direction  (figure 
4.8)  shows  excellent  agreement  between  the  continuum  model  and  VCFEM-HOMO  results. 
An  increase  in  the  inclusion  volume  fraction  leads  to  a  considerable  reduction  in  the  ductility 
due  to  the  higher  stiffness  of  the  elastic  inclusion.  Though  the  random  distribution  shows 
a  less  stiff  response  than  the  square  edge  distribution  for  the  same  volume  fraction  (55%), 
this  plot  confirms  that  the  effect  of  volume  fraction  is  much  more  pronounced  that  that  of 
distribution.  The  two  approaches  are  at  considerable  variance  in  the  transverse-longitudinal 
strain  plot  (figure  4.9).  The  continuum  model  shows  nearly  identical  strain  responses  for 
all  inclusion  volume  fractions  and  distributions.  This  can  be  explained  by  the  fact  that  the 
continuum  models  assume  total  yielding  of  the  microstructure  after  initial  yield,  and  conse¬ 
quently  enforces  the  constraint  that  macroscopic  plastic  strains  are  volume  preserving.  This 
leads  to  direct  dependence  of  transverse  plastic  strains  on  the  longitudinal  strains,  and  the 
microstructure  contributes  only  to  the  elastic  part  of  the  transverse  strain.  The  behavior 
is  however  quite  different  when  solved  with  VCFEM-HOMO.  In  reality,  only  a  part  of  the 
microstructure  yields  at  the  strains  considered,  and  this  is  reflected  in  the  VCFEM-HOMO 
results.  For  smaller  inclusion  volume  fractions,  a  large  portion  of  the  unit  cell  is  the  matrix 
that  goes  plastic  at  initial  yield.  This  gives  the  close  proximity  in  prediction  by  both  models. 
This  behavior  changes  significantly  at  higher  volume  fraction,  where  larger  portions  of  the 
matrix  remain  elastic.  The  absolute  value  of  the  transverse  strain  by  VCFEM-HOMO  is 
much  smaller  than  that  of  effective  continuum  model.  Another  observation  is  that  the  differ¬ 
ence  between  VCFEM-HOMO  and  continuum  model  results  for  the  random  microstructure  is 
less  than  that  for  the  square  edge  microstructure.  This  may  be  attributed  to  the  anisotropic 
behavior  of  the  latter  that  is  not  accounted  for  in  the  continuum  models.  This  example 
clearly  points  out  the  need  for  using  the  coupled  analysis  especially  at  higher  volume  frac¬ 
tions  of  heterogeneous  microstructures. 

The  second  example  deals  with  porous  materials,  for  which  VCFEM-HOMO  results  are 
compared  with  the  Tvergaard-Gurson  (T-G)  continuum  model  [38,  33].  Plane  strain  simple 
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Figure  4.8:  Macroscopic  stress-strain  relations  for  various  microstructures  by  VCFEM-HOMO  and 
effective  continuum  model 
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Macroscopic  strain  in  loading  direction 


Figure  4.9:  Transverse  strain-axial  strain  relations  for  various  microstructures  by  VCFEM-HOMO 
and  effective  continuum  model 
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tension,  simple  shear,  and  biaxial  tension  tests  are  carried  out  for  square  edge  packed  ma¬ 
terials  with  a  low  void  volume  fraction  of  5%.  For  comparison  plane  strain  simulations,  the 
T-G  model  for  long  cylindrical  voids  [38]  is  implemented  in  ABAQUS  through  the  UMAX 
window.  In  this  module,  numerical  integration  of  the  rates  of  state  variables  is  carried  out  in 
a  way  similar  to  Aravas  [40],  for  spherical  voids.  The  macroscopic  ABAQUS  model  for  both 
analyses  consists  of  one  QUAD4  element  and  the  microscopic  model  for  VCFEM-HOMO  is 
illustrated  in  figure  4.3.  The  yield  condition  in  the  plane  strain  T-G  model  [38]  is  stated  as  : 


$  = 


a 


+  3  /  cosh{- 


\/3  (<Tii  -f-  <722) . 


-  (l+2.25f) 


(4.52) 


where  is  the  effective  stress,  do  is  the  matrix  flow  stress,  /  is  the  void  volume  fraction  and 
dij  are  the  principal  components  of  the  Cauchy  stress  tensor.  An  important  assumption  in 
this  model  is: 

(1  -  f)<^od  =  (TqeU  (4.53) 

which  means  that  the  equivalent  plastic  work  is  derived  solely  from  the  entire  matrix  part 
of  the  microstructure. 


Effective  stress-strain  plots  in  simple  tension  and  simple  shear  (figures  4.10(a)  and  (b)) 
show  good  agreement  between  two  methods,  with  T-G  model  producing  a  slightly  higher 
value  of  initial  yield.  This  is  because,  yield  in  the  T-G  model  is  manifested  only  after  a 
significant  portion  of  the  microstructure  has  become  plastic.  VCEEM-HOMO  on  the  other 
hand  shows  signs  of  inelasticity  at  the  very  onset  of  plastic  deformation  in  the  microstructure. 
The  two  approaches  are  however  at  significant  variance  for  the  biaxial  loading  case  as  seen 
in  figure  4.11(a).  Though  the  initial  yield  points  are  fairly  close,  the  post  yield  behavior 
shows  considerable  difference  between  two  models.  As  in  the  previous  example,  the  main 
reason  for  this  difference  may  be  attributed  to  the  assumption  in  T-G  model  that  the  entire 
microstructural  matrix  becomes  plastic  upon  initial  yield.  In  biaxial  loading,  as  {dn  -h  (T22) 
increases  in  equation  (4.52),  the  effective  macroscropic  stress  cr^  must  diminish  to  near  zero 
values  for  satisfying  the  yield  condition.  Physically,  this  also  indicates  that  the  large  regions 
of  plastic  localization  causes  the  effective  stress  to  drop  drastically.  Erom  the  VCEEM- 
HOMO  simulation  it  is  however  seen  in  figure  4.11(b),  that  a  significant  portion  of  the 
microstructure  does  not  yield.  Consequently,  the  effective  stress  in  figure  4.11(a)  is  higher 
by  VCEEM-HOMO  than  by  the  T-G  predictions  at  the  end  of  loading. 

As  a  final  example  a  multiple  scale  analysis  for  a  porous  plate  with  a  hole  is  conducted 
by  VCEEM-HOMO  and  compared  with  the  T-G  model  incorporated  in  ABAQUS.  Uniaxial 
tension  of  the  plate  is  considered  under  plane  strain  conditions.  The  representative  material 
element  (RME)  in  the  VCFEM-HOMO  analysis  is  a  square  edge  packing  with  5%  circular 
void  volume  fraction.  Eigures  4.12(a)  and  (b)  shows  excellent  agreement  in  contour  plots  of 
effective  macroscopic  stress  distribution  by  VCFEM-HOMO  and  ABAQUS  with  T-G  model. 
The  stress-strain  evolution  at  three  specific  points  (A),  (B)  and  (C)  are  compared  in  figure 
4.12(c).  This  near  perfect  match  is  expected,  since  most  points  in  the  plate  are  essentially 
in  an  uniaxial  state  of  stress. 
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Macroscopic  effective  stress(Gpj 


Figure  4.10:  Effective  macroscopic  stress-strain  plots  by  VCFEM-HOMO  and  Tvergaard-Gurson 
model  for  5%  void  volume  fraction  in  plane  strain  (a)  simple  tension  and  (b)  simple  shear 


Macroscopic  effective  stress(GPa) 


(a) 


Figure  4.11:  (a)  Effective  macroscopic  stress-strain  plots  by  VCFEM-HOMO  and  Tvergaard-Gurson 
model  for  5%  void  volume  fraction  in  biaxial  tension  and  (b)  Contour  plot  of  microstructural 
effective  plastic  strain  at  the  end  of  loading 
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Comparisons  with  some  analytical  and  experimental  results 

Plane  strain  VCFEM-HOMO  results  are  compared  with  analytical  2  and  3-phase  model 
predictions  due  to  Zhao  and  Weng  [41]  and  experimental  results  due  to  Adams  [42]  for 
unidirectional  composites.  The  material  properties  are: 

Aluminum 

E  =  58  GPa  (Young's  Modulus  ),  i/  =  0.33  (Poisson's  Ratio) 

Yo  =  89 -¥Pa  (Initial  Yield  Stress),  (Tegt;  =  io  + 175(eP^„)°-^2^  MPa  (Post  Yield  Hardening 
Law  ) 

Boron 

E  =  385  GPA  (Young’s  Modulus  ),  u  —  0.2  (Poisson’s  Ratio) 

A  square  macroscopic  element,  that  has  a  microstructural  RME  of  a  2  x  1  rectangular  edge 
packing  with  a  circular  inclusion  of  Vj  =  34%,  is  loaded  in  uniaxial  tension.  Figure  4.13 
shows  that  VCFEM-HOMO  results  provide  a  better  agreement  with  experimental  results 
than  the  analytical  model,  which  assume  a  square  edge  packed  RME. 

Characteristic  modes  in  elastic-plastic  deformation 

Characteristic  modes  in  elastic-plastic  homogenization  represent  instantaneous  microstruc¬ 
tural  deformation  response  to  applied  unit  components  of  macroscopic  strain  tensor.  These 
modes  Xp^  are  obtained  from  equation  (4.49).  Characteristic  modes  are  essential  in  obtaining 
the  elastic-plastic  tangent  modulus  EPf.i  for  macroscopic  analysis. 

The  macroscopic  element  is  subjected  to  uniaxial  stretching,  and  the  microstructural 
RMEs  considered  are  square  edge/random  packing  for  circular  voids  (20%  volume  fraction), 
and  square  edge  packing  for  circular  inclusion  (40%  volume  fraction).  Figure  4.14(a)  and  (b) 
show  the  characteristic  modes  for  the  porous  material  and  corresponding  tangent  modulus 
at  the  beginning  of  the  loading  and  after  1%  macroscopic  stretching  respectively.  The  initial 
modes  depict  the  deformed  void  aligned  in  the  direction  of  straining,  but  this  changes  with 
plastic  straining,  implying  the  influence  of  internal  variable  evolution.  Figure  4.15  shows 
corresponding  modes  for  the  random  microstructure.  Figure  4.16  shows  characteristic  modes 
for  composite  microstructure  with  regular  distribution.  The  difference  between  the  modes 
for  porous  and  composite  materials  emanates  from  the  high  inclusion  stiffness  that  results 
in  larger  matrix  deformation. 

4.6.2  Multiple  Scale  Analysis 

Two  problems  are  considered  for  multiple  scale  analysis  with  VCFEM-HOMO.  The  effects 
of  various  distributions,  sizes,  and  shapes,  on  the  overall  behavior  of  the  heterogeneous 
material  is  investigated.  In  particular,  the  RME’s  considered  are  (a)  square  edge  packing 
with  a  circular  heterogeneity,  (b)  square  diagonal  packing  with  a  circular  heterogeneity,  (c) 
random  packing  with  15  circular  heterogeneities  and  (d)  random  packing  with  15  elliptical 
heterogeneities.  The  material  microstructure  is  assumed  to  consist  of  voids,  inclusions,  or  a 
combination  of  both. 
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Macroscopic  tensile  plastic  strain{%) 


Figure  4.13:  Macroscopic  stress-strain  responses  for  a  Boron- Aluminum  composite  by  VCFEM- 
HOMO,  analytical  and  experimental  techniques 
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Figure  4.14;  Characteristic  modes  for  voided  RME  with  square  edge  distribution,  (a)  at  the  start 
of  loading  (b)  after  1%  stretching 
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Figure  4.15:  Characteristic  modes  for  voided  RME  with  random  distribution,  (a)  at  the  start  of 
loading  (b)  after  1%  stretching 
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Figure  4.16;  Characteristic  inodes  for  composite  RME  with  square  edge  distribution,  (a)  at  the 
start  of  loading  (b)  after  1%  stretching 


Plane  strain  analysis  of  plate  with  holes 

This  set  of  examples  deals  with  a  thick  plate  with  an  uniform  array  of  large  circular  holes  in 
plane  strain  uniaxial  tension.  Only  a  portion  of  the  plate,  shown  in  figure  4.17(a),  is  consid¬ 
ered  from  symmetry  considerations.  Dimensional  details  are  given  in  the  figure  4.17(a).  The 
macroscopic  ABAQUS  model  consists  of  128  QUAD4  elements.  The  left  and  right  edges  are 
constrained  to  move  in  vertical  straight  lines,  and  the  top  and  bottom  edges  are  pulled  to 
an  overall  strain  of  0.5%. 

In  the  first  example,  a  porous  microstructure  with  a  40%  void  volume  fraction  is  consid¬ 
ered.  Various  morphologies  are  considered  are  shown  in  figures  4.17  (b),  (c),  (d)  and  (e). 
Figures  4.18,  4.20,  4.22  and  4.24  show  the  effective  stress  contours  plots,  and  figures  4.19, 
4.21,  4.23  and  4.25  show  the  contours  plots  of  the  effective  plastic  strains  at  the  structural 
level  and  also  in  a  microstructural  RME  at  the  corner  point  A.  It  is  observed  that  the 
maximum  effective  stress  generally  occurs  at  this  point.  The  macroscopic  effective  stress 
distributions  exhibit  very  similar  patterns  for  all  microstructures  considered.  There  is  a 
narrow  ligament  between  two  large  holes  along  which  strains  localize  as  the  deformation 
intensifies.  Similarity  in  the  patterns  of  macroscopic  stress  distribution  is  explained  from  an 
observation  that  the  macroscopic  state  of  stress  at  each  point  is  essentially  uniaxial.  Thus 
anisotropy  emanating  from  the  microstructural  morphology  does  not  play  an  important  role 
in  this  example.  Magnitude  of  stresses  and  plastic  strains  are  however  considerably  different 
depending  on  the  microstructural  arrangement.  This  is  further  evidenced  in  a  macroscopic 
stress-strain  plot  (figure  4.26  (a))  at  the  point  A  .  The  square  edge  and  square  diagonal 
distributions  yield  nearly  identical  response.  The  random  void  distributions  exhibit  signif¬ 
icantly  more  ductile  behavior  compared  to  the  regular  distributions.  The  ductility  in  the 
random  distribution  stems  from  the  plastic  strain  concentration  in  regions  of  clustered  voids. 

A  comparison  of  the  macroscopic  and  microscopic  contour  plots  reveals  that  the  true 
stress  in  the  microstructure  is  significantly  higher  than  the  macroscopic  stresses.  For  exam¬ 
ple,  at  point  A  the  maximum  microscopic  effective  stress  is  160%  higher  for  square  edge, 
175%  higher  for  square  diagonal,  265%  higher  for  random  circular  and  340%  higher  for  ran¬ 
dom  elliptical  packing.  An  interesting  observation  is  that,  though  the  random  distribution 
models  shows  more  ductile  behavior,  it  gives  rise  to  local  regions  of  significantly  higher  mi¬ 
croscopic  stresses.  In  the  analysis  of  porous  materials  with  large  void  volume  fractions,  the 
interaction  between  microstructural  voids,  and  consequently  the  inter- void  distances  appears 
to  play  a  more  dominant  role  than  the  arrangement  itself,  both  for  macroscopic  and  micro¬ 
scopic  response. 

The  same  problem  is  considered  again  with  the  microstructure  changed  to  a  boron- 
aluminum  composite.  The  same  microstructural  morphologies  are  used  with  40%  fiber 
volume  fraction.  Figures  4.27,  4.28,  4.29  and  4.30  show  the  effective  stress  contour  plots 
at  the  end  of  loading.  The  macroscopic  plots  show  that  the  localization  ligament  between 
two  macroscopic  holes  is  much  less  severe  and  more  diffused  compared  to  the  porous  ma¬ 
terial.  For  the  same  level  of  macroscopic  strain,  the  square  edge  packing  yields  the  highest 
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values  of  effective  stress  while  the  square  diagonal  packing  yields  the  lowest,  with  the  ran¬ 
dom  packings  in  between.  This  is  also  evidenced  in  the  macroscopic  stress-strain  plot  at 
point  A  in  figure  4.26(b).  As  opposed  to  the  porous  microstructure,  the  arrangement  of  the 
composite  microstructure  has  a  major  effect  on  the  macroscopic  response.  The  distribution 
pattern  is  again  similar  for  all  RME’s  due  to  the  predominantly  uniaxial  state  of  local  stress. 
The  microscopic  contour  plots  reveal  that  the  true  stress  in  the  microstructure  is  signifi¬ 
cantly  higher  than  the  macroscopic  counterparts.  For  example,  at  point  A  the  maximum 
microscopic  effective  stress  is  91%  higher  for  square  edge,  50%  higher  for  square  diagonal, 
180%  higher  for  random  circular  and  177%  higher  for  random  elliptical  packing.  This  shows 
that  high  stress  concentration  occurs  in  certain  fibers  for  random  microstructures,  and  are 
consequently  more  susceptible  to  microstructural  damage  by  fiber  cracking  or  debonding  for 
identical  levels  of  macroscopic  stresses. 


Thick  cylinder  with  internal  pressure 

In  this  example  a  thick  composite  cylinder  is  subjected  to  internal  pressure  that  is  increased 
from  0  to  a  maximum  of  0.1  GPa.  Only  a  quarter  of  the  cylinder  is  modeled  from  symme¬ 
try  considerations  with  60  QUAD4  elements  in  ABAQUS.  The  dimensions  of  the  cylinder 
and  boundary  conditions  are  shown  in  figure  4.31(a).  The  cylinder  is  analyzed  for  three 
different  microstructures  with  40%  fiber  volume  fraction  under  plane  strain  conditions.  The 
microstructural  VCFEM  models  for  square  edge,  square  diagonal  and  random  packings  are 
displayed  in  figures  4.31(b),  (c)  and  (d). 

Figures  4.32,  4.34,  4.36  show  contour  plots  of  the  Von  Mises  stress  and  4.33,  4.35,  4.37 
show  the  corresponding  effective  plastic  strains  at  macroscopic  and  microscopic  levels  at 
the  end  of  loading.  The  macroscopic  plots  show  that  for  the  random  packing  the  stresses 
and  strains  depend  only  on  the  radial  location  (r)  and  not  the  angular  location  {6)  of  the 
point.  However,  for  the  square  edge  and  square  diagonal  packings,  these  state  variables 
depend  on  both  the  (r)  and  [0)  coordinates,  leading  to  "earing”  effects.  The  results  of 
the  square  edge  and  square  diagonal  packings  are  mutually  rotated  by  45°  as  expected 
from  the  microstructural  orientations.  Such  earing  effects  are  known  to  be  consequences 
of  anisotropy  in  the  microstructures.  In  this  problem,  most  macroscopic  points  are  in  a 
triaxial  state  of  state,  as  opposed  to  the  previous  example  where  most  points  were  in  a 
uniaxial  stress  state.  This  amplifies  the  effect  of  anisotropy  in  the  macroscopic  stress  and 
strain  distributions.  Effective  behavior  of  the  random  packing  is  near  isotropic  and  hence 
the  stress  strain  distributions  are  independent  of  9  coordinates. 
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Figure  4.18:  Von  Mises  stress  distribution  for  porous  material  with  square  edge  packing  (a)  macro- 
scopic  stress,  (b)  microscopic  stress  at  point  A 
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Figure  4.23:  Effective  plastic  strain  distribution  for  porous  material  with  random  circular  packing 
(a)  macroscopic  strain,  (b)  microscopic  strain  at  point  A 
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Figure  4.25:  Effective  plastic  strain  distribution  for  porous  material  with  random  elliptical  packing 
(a)  macroscopic  strain,  (b)  microscopic  strain  at  point  A 
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Figure  4.26:  Evolution  of  effective  stress  with  evg^ving  strain  at  point  A  in  the  heterogeneous  plate 
with  40%  second  phase  volume  fraction  for  (a)  porous  material  (b)  composite  material 
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F'igure  4.33:  Effective  plastic  strain  distribution  for  square  edge  packing  (a)  macroscopic  stress,  (b) 
microscopic  stress  at  point  A 


4*7  Conclusions 


In  this  work,  a  multiple-scale  computational  tool  (VCFEM-HOMO)  is  devised  for  performing 
elastic-plastic  analysis  of  heterogeneous  materials  with  inclusions  and  voids  in  the  microstruc¬ 
ture.  The  microscopic  analysis  is  conducted  with  the  Voronoi  Cell  finite  element  model  while 
a  conventional  displacement  based  FEM  code  (ABAQUS  in  this  work)  executes  the  macro¬ 
scopic  analysis.  Coupling  between  the  scales  is  accomplished  through  the  user  based  UMAT 
window  in  ABy\QUS,  for  incorporating  the  effective  constitutive  model  from  VCFEM  by 
asymptotic  homogenization. 

The  VCFEM  developed  for  porous  and  composite  materials  in  [26,  27]  is  used  in  conjunc¬ 
tion  with  homogenization  for  a  wide  variety  of  arbitrary  microscopic  phase  distributions. 
In  this  model  for  small  deformation  elasto-plasticity,  stress  functions  are  motivated  from 
essential  characteristics  of  micromechanics  by  accounting  for  the  heterogeneity  shape  and 
influence.  This  is  achieved  by  transforming  any  arbitrary  shaped  heterogeneity  to  a  unit  cir¬ 
cle  using  a  mapping  similar  to  the  Schwarz-Christoffel  conformal  mapping.  Stress  functions 
are  then  constructed  in  terms  of  mapped  coordinates.  The  introduction  of  micromechanics 
consideration  tremendously  enhances  VCFEM  accuracy  for  various  microstructures  at  very 
moderate  computational  efforts.  The  accuracy  and  efficiency  of  VCFEM  are  established 
by  comparing  with  conventional  FEM  commercial  packages.  For  a  wide  range  of  problems 
VCFEM  delivers  very  similar  accuracy  at  a  considerably  low  computational  effort.  This  is 
evidenced  by  the  drastically  reduced  degrees  of  freedom  needed  for  convergence,  compared 
with  the  conventional  codes.  The  D.O.F.  ratio  varies  from  as  low  as  20  for  simple  microstruc¬ 
tures.  to  even  100  for  more  complex  cases.  This  translates  into  a  reduction  factor  of  15-30 
in  the  CPU  time  for  execution,  even  with  a  non-optimized  research  code.  Furthermore,  user 
effort  required  to  generate  the  model  is  much  lower  for  VCF"EM  than  for  many  commercial 
packages. 

Numerical  examples  conducted  with  VCFEM-HOMO  establish  the  effectiveness  of  ho¬ 
mogenization,  when  compared  with  the  FEM  calculations  with  constitutive  relations  from 
unit  cell  and  effective  continuum  models.  They  also  point  out  the  limitations  of  assumptions 
made  in  the  latter  methods,  and  emphasize  situations  when  coupled  multiple  scale  analy¬ 
sis  becomes  necessary.  The  effect  of  various  microscopic  arrangements  on  the  mechanical 
response  at  the  two  scales  are  investigated  through  these  examples.  Significant  influence  is 
observed  in  both  the  microscopic  and  macroscopic  responses.  For  example,  the  square  edge 
and  square  diagonal  distributions  give  severe  anisotropic  effects,  while  random  distributions 
are  more  isotropic  at  a  point  in  the  structure.  Multiple  scale  computations  with  micro¬ 
scopic  FEM  models  at  each  element  integration  point  in  the  structure  are  in  general  quite 
exhaustive.  For  complex  microstructures,  the  efficiency  of  microstructural  VCFEM  makes 
it  possible  to  realize  computations  of  this  magnitude  in  comparison  with  conventional  FEM. 
The  more  complex  multiple  scale  models  in  this  work  took  approximately  2-3  CPU  hours 
to  run  on  a  CRAY-YMP  at  the  Ohio  Supercomputer  Center.  Though  reasonable,  major 
improvements  in  efficiency  of  the  scale  coupling  is  presently  being  pursued. 
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In  conclusion,  the  Voronoi  Cell  finite  element  model  with  asymptotic  homogenization 
(VCFEM-HOMO)  emerges  as  an  important  tool  for  analyzing  arbitrary  microstructures  in 
many  materials.  It  is  easily  adapted  with  commercial  packages  at  the  structural  scale,  and 
this  makes  it  very  attractive.  A  shortcoming  of  the  homogenization  method  in  its  present 
form,  is  however  the  incorporation  of  boundary  condition  at  regions  of  material  discontinuity. 
Boundary  conditions  should  be  applied  in  the  microstructure  and  not  at  the  macroscopic 
level  as  is  presently  done.  The  boundary  effect  can  become  pronounced  in  some  cases,  and 
then  homogenization  results  become  less  accurate  near  the  boundary  (see  [22,  23]). 
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Chapter  5 


R-S  Adapted  Arbitrary 
Lagrangian-Eulerian  Finite  Element 
Method  for  Metal  Forming  Problems 
with  Strain  Localization 

Summary 

In  this  chapter,  an  adaptive  arbitrary  Lagrangian-Eulerian  (ALE)  finite  element  method 
is  developed  for  solving  large  deformation  problems,  with  applications  in  metal  forming 
simulation.  The  ALE  mesh  movement  is  coupled  with  r-adaptation  of  automatic  node  relo¬ 
cation,  to  minimize  element  distortion  during  the  process  of  deformation.  Strain  localization 
is  considered  in  this  study  through  the  constitutive  relations  for  ductile  porous  materials, 
proposed  by  Gurson  and  Tvergaard.  Prediction  of  localized  deformation  is  achieved  through 
a  multi-level  mesh  superimposition  method,  termed  as  s-adaptation.  The  model  is  validated 
by  comparison  with  established  results  and  codes,  and  a  few  metal  forming  problems  are 
simulated  to  test  its  effectiveness. 


5.1  Introduction 

Metal  forming  processes  involve  large  plastic  deformation  of  a  workpiece  to  conform  to  the 
desired  shape  of  the  final  product.  For  certain  materials  and  processes  severe  localized  defor¬ 
mation  can  lead  to  the  formation  of  localization  bands  that  may  render  a  product  defective. 
It  is  very  important  to  accurately  simulate  such  large  strain  problems  and  redesign  the  pro¬ 
cess  to  avoid  the  existence  of  these  defects  in  the  product.  This  paper  is  aimed  at  devising 
an  accurate  and  efficient  adaptive  finite  element  model  to  simulate  metal  forming  processes 
where  there  is  a  possibility  of  localized  straining. 

Such  a  simulation  offers  many  challenges  which  should  be  systematically  overcome  to 
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minimize  user  intervention  and  experience.  For  example,  these  simulations  often  lead  to 
severe  mesh  distortions  and  inadequate  representation  of  the  contact  boundary  for  sharp 
edges  in  the  worktools.  The  arbitrary  Lagrangian-Eulerian  (ALE)  kinematic  description  has 
been  successfully  implemented  in  analysis  of  large  deformation  in  solids  [1,  2,  3,  4]  to  avoid 
these  difficulties.  In  this  description,  the  finite  element  mesh  is  decoupled  from  the  under¬ 
lying  material,  thereby  creating  a  flexible  mesh  which  can  have  a  motion  different  from  the 
material.  The  ALE  description  is  used  to  overcome  some  of  the  problems  encountered  with 
a  pure  Lagrangian  or  a  pure  Eulerian  type  of  mesh.  Applications  of  the  ALE  description  in 
metal  forming  problems  have  been  examined  by  Ghosh  et.al.  [3,  4,  5],  Huetink  et.  al.  [2], 
Liu  et.  al.  [6]  among  others. 

Despite  the  flexibility  of  the  finite  element  mesh  in  the  ALE  description,  additional  de¬ 
grees  of  freedom  introduced  can  pose  difficulties  while  simulating  complex  problems.  Con¬ 
trolling  the  motion  of  an  ALE  mesh  requires  a  keen  knowledge  on  the  deformation  process, 
and  can  adversely  affect  the  simulation  if  improperly  specified.  Lack  of  knowledge  on  the 
sensitivity  of  solutions  to  the  grid  configuration  can  deteriorate  the  quality  of  solutions.  Ob¬ 
stacles  in  the  ALE  method  can  be  surpassed  by  combining  with  grid  adaptation  methods  as 
a  tool  for  reducing  element  distortions  incurred  through  improper  nodal  positioning  during 
deformation.  The  r-method  of  node  relocation  driven  by  element  distortion  is  consequently 
coupled  with  the  ALE  description.  This  method  has  been  applied  to  large  deformation  prob¬ 
lems  by  Kikuchi  and  Torigaki  [7]  to  reduce  interpolation  error. 

Strain  localization  occur  over  very  narrow  regions  in  comparison  with  the  overall  di¬ 
mensions  of  the  workpiece.  Since  most  constitutive  relations  are  local  in  nature,  the  finite 
element  predictions  are  extremely  mesh  sensitive.  Consequently  for  many  meshes  with  no 
apriori  knowledge  of  localized  regions,  a  simulation  may  not  be  capable  to  capture  this  phe¬ 
nomenon  due  to  lack  of  appropriate  resolution.  A  coarse  mesh  may  completely  miss  the 
localized  behavior  in  deformation.  Various  computational  methods  have  evolved  for  simu¬ 
lating  problems  of  strain  localization  without  having  to  use  extremely  high  resolution  in  the 
overall  mesh.  Important  contributions  among  these  include  coupled  finite  element-spectral 
methods  by  Belytschko  et.  al.  [23,  24],  assumed  stress  based  finite  element  by  Atluri  et  al 
[25],  [26],  finite  element  with  embedded  localization  zones  and  the  enhanced  strain  elements 
by  Simo  and  Armero  [27,  28]. 

Adaptive  finite  element  techniques  have  evolved  to  optimize  the  overall  computational 
effort  and  provide  a  measure  of  reliability  through  error  indicators.  An  important  feature  of 
the  adaptive  method  is  to  use  computed  results  a-posteriori  for  improving  the  quality  of  the 
solution.  Various  methods  such  as  h-method  of  mesh  refinement  [8,  11,  12,  14],  and  p-method 
of  mesh  enrichment  [15,  16]  and  r-method  of  node  relocation  [17]  have  been  proposed.  Of 
all  methods  the  combined  h-p  method  adaptation  have  been  proven  to  be  very  effective  for 
a  wide  class  of  problems  with  localization  such  as  shock  wave  propagation  and  strain  local¬ 
ization  etc  [20,  19,  13].  Fish  et.  al.  [21,  22],  have  proposed  a  powerful  adaptive  multi-level 
mesh  superposition  method  called  the  s-method.  This  method  increases  the  resolution  of  the 
finite  element  analysis  by  superimposition  of  multiple  levels  of  hierarchical  elements  in  local 
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regions  of  high  solution  gradients.  Fish  and  his  coworkers  have  applied  this  method  for  a 
wide  variety  of  problems,  including  linear  elastic  analysis  of  homogeneous  and  heterogeneous 
materials,  laminated  composites,  and  problems  with  cracks. 

In  this  paper,  the  r-adapted  ALE  mesh  is  combined  with  multilevel  s-adaptation  to  sim¬ 
ulate  metal  forming  problems  featuring  large  elastic-plastic  deformation  with  strain  localiza¬ 
tion.  Though  the  global  mesh  is  ALE,  the  s-adapted  mesh  is  made  to  follow  the  material.  The 
material  yield  function  is  that  for  porous  ductile  materials,  originally  proposed  by  Gurson 
[31]  and  subsequently  modified  by  Tvergaard  [34,  36].  A  perturbed  Lagrangian  formulation 
is  executed  for  normal  contact  between  workpiece  and  worktools. 


5.2  Governing  Equations 

5.2.1  Basic  equations  in  the  ALE  model 

The  arbitrary  Lagrangian-Eulerian  description  introduces  a  reference  configuration,  consist¬ 
ing  of  a  set  of  unconstrained  grid  points  that  undergo  arbitrary  spatial  motion.  Each  point 
in  this  reference  configuration  is  unambiguously  represented  by  an  invariant  set  of  three  in¬ 
dependent  coordinates  Xi-  The  motion  of  points  in  the  reference  frame  is  denoted  by  the  set 
of  spatial  coordinates  x,,  and  is  expressed  as  arbitrary  continuous  functions  of  Xi  and  time 
t.  The  continuity  and  momentum  equations  in  this  description  are  respectively  written  as: 
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where  W  refers  to  the  velocity  of  a  representative  grid  point,  V  is  the  velocity  of  the  material 
point  coinciding  with  the  grid  at  time  t,  cTij  is  the  Cauchy  stress  tensor,  p  is  the  density,  6, 
is  the  body  force  per  unit  volume  and  J  is  the  Jacobian.  Evaluation  of  process  variables 
at  nodal  points  in  an  ALE  description  requires  an  additional  relation  for  updating  variable 
values  from  material  points  to  grid  points.  This  equation  correlates  time  derivatives  of  a 
material  variable  in  fixed  material  coordinates  to  that  in  fixed  referential  coordinates.  Thus, 
if  ^  is  a  time  dependent  material  variable,  rates  in  the  two  coordinate  systems  are  related 

3<S  * 
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where  [^  and  ]^  refer  to  material  and  referential  time  derivatives  respectively. 


5.2.2  Constitutive  Model 

Material  parameters  play  a  very  important  role  in  the  development  of  localized  shear  bands 
when  subjected  to  quasi-static  loading  conditions.  Localization  has  been  proven  to  occur 
with  features  such  as  yield  surface  curvature  [34],  vertex  elfects  [35],  alignment  of  preferred 
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slip  direction  [40]  and  internal  damage  accumulation  such  as  void  volume  fraction  [34]  etc. 
Since  the  metal  forming  simulations  considered  in  this  study  are  conducted  under  quasi¬ 
static  conditions,  a  constitutive  model  featuring  void  nucleation,  growth  and  coalescence  is 
adopted.  This  model  has  been  developed  for  porous  ductile  materials  by  Gurson  [31]  for 
void  growth,  and  subsequently  modified  by  researchers  like  Tvergarrd  [34,  36],  Tvergaard 
and  Needleman  [35]  etc  for  void  growth  and  coalescence.  Void  nucleation  model  in  this  re¬ 
search  is  based  on  the  model  of  Needleman  and  Rice  [39].  Horn  and  McMeeking  [37]  have 
shown  that  the  Gurson-Tvergaard  model  compares  very  well  with  complete  3-D  finite  ele¬ 
ment  predictions  for  uniaxial  tension,  hydrostatic  tension,  and  pure  shear  stress  fields. 

In  this  model  an  approximate  yield  condition  is  postulated  in  the  form 
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where  is  the  macroscopic  effective  stress,  crjis  the  equivalent  tensile  flow  stress  representing 
actual  microscopic  stress  state  in  the  matrix  material,  (Tkk  is  the  hydrostatic  stress  and  /  is 
the  void  volume  fraction.  The  control  function  /*  is  defined  in  terms  of  /  as 
/*(/)  =  /  for  f<fc 

=  fc  +  j^ifF  -  fc)  for  f>fc 
where  fc  is  a  critical  value  of  void  volume  fraction  (~  0.15),  fu  is  the  ultimate  value  of  void 
volume  fraction  at  fracture  (~  0.25),  and  =  /(/f)  =  ^  is  the  volume  fraction  at  which 
the  material  loses  it  stress  carrying  capacity  [39].  The  parameters  ^i,  q2  and  qs  are  chosen  to 
be  1.5,  1.0  and  2.25  respectively  in  this  study  as  suggested  in  [33,  36].  The  rate  of  evolution 
of  /  is  taken  to  be 
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Void  nucleation  parameters  A  and  B,  can  be  obtained  for  cases  like  plastic  strain  controlled, 
stress  controlled  or  both,  as  suggested  in  [39].  In  this  work  only  plastic  strain  controlled 
nucleation  is  considered  for  which 


E,  Es^/Wi 


P  2 


B  =  ^ 


In  equation  (5.6),  Et  corresponds  to  the  elastic-plastic  tangent  modulus,  f^  is  the  volume 
fraction  of  void  nucleating  particles,  cat  is  the  mean  nucleation  strain  and  s  is  the  correspond¬ 
ing  standard  deviation.  Based  on  [32],  values  /at  =  0.04,  s  =  0.1  and  =  0.3  are  used  in 
this  study.  The  plastic  part  of  the  macroscopic  rate  of  deformation  tensor  is  obtained 
by  normality  rule  as 


An  additional  assumption  that  the  rate  of  equivalent  work  in  matrix  material  equals  the 
macroscopic  work  is  made  as: 
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=  (!-/)  <’1  i)  (5-8) 

Finally,  the  rate  independent  elastic-plastic  equation  is  formulated  in  a  rotated  Lagrangian 
system  from  considerations  of  stress-rate  objectivity  as  explained  in  [3,  4]  and  is  expressed 
as 


tij  =  Eijki{dki  —  d^i)  (5.9) 

where  Uj  is  material  derivative  of  a  rotated  Cauchy  stress  tensor  Uj,  obtained  by  rotation 
of  the  Cauchy  stress  tensor  aij  through  Rij  which  is  derived  from  the  polar  decomposition 
of  the  deformation  gradient  tensor,  dij  is  the  corresponding  rotated  rate  of  deformation 
tensor,  and  Eijki  is  a  fourth  order  elasticity  tensor.  Numerical  integration  of  the  elasto- 
plastic  constitutive  equations  described  above,  is  carried  out  using  a  technique  based  on  the 
backward  Euler  integration  rule,  proposed  by  Aravas  [38]. 

5.2.3  Contact  Formulation 

Only  normal  contact  is  considered  in  this  study.  The  perturbed  Lagrangian  formulation  ( 
PLF  )  [41]  is  implemented  to  incorporate  non-penetration  constraint  on  the  contact  surface. 
This  method  combines  the  advantages  of  Lagrange  multiplier  method  and  the  penalty  func¬ 
tion  method  to  improve  convergence  and  reduce  the  dependence  on  penalty  parameter.  It 
introduces  a  regularization  of  the  Lagrange  multiplier  terms  by  a  penalty  term  to  a  modified 
form  of  the  potential  energy  as: 

n,  =  npe+/  X<g>dT-^l  A^dF  (5.10) 

where  Ilpe  is  the  strain  energy  of  the  deforming  body,  A  is  the  Lagrange  multiplier,  ^  is  a 
gap  function  representing  distance  between  contacting  surfaces  and  F*^  is  the  contact  surface, 
e  is  an  arbitrary  positive  parameter  and  <>  is  the  MacCaulay  operator.  Equation  (  5.10) 
approaches  the  Lagrange  multiplier  method  as  e  tends  to  infinity.  The  contact  pressure  is 
assumed  to  be  constant  in  each  contact  segment,  yielding 

Xs  =  -{9i+92)  (5.11) 

Here  ^ri  and  92  are  the  gaps  at  the  edges  of  the  contact  segment  as  shown  in  figure  5.1. 
This  results  in  an  averaged  form  of  the  non-penetration  constraint  (  >  0)  for  each  contact 

segment. 

During  the  deformation  process,  when  two  surfaces  come  in  contact  as  shown  in  figure  5.1, 
the  following  steps  are  implemented  to  incorporate  the  contact  condition  (5.11). 

1.  Drop  a  perpendicular  from  node  (B)  of  surface  1  onto  segment  (ab)  of  surface  2  and 
identify  the  intersection  (f). 

2.  Similarly  drop  a  perpendicular  from  node  (b)  on  surface  2  onto  segment  (BC)  of  surface 
1  and  identify  the  intersection  (F). 
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Surface  1 


Surface  2 
e 


'  contact  segnients^ 

Figure  5.1:  Discretization  of  the  contact  interface  into  contact  segments. 

Contact  segments  for  surface  2  are  fb,  bg,  gc,  ch,  etc.  The  corresponding  contact  segments 
for  surface  1  are  BF,  FC,  CG,  GD  etc.  The  gap  functions  gi,  ...  are  shown  in  figure  5.1 
from  which  contact  pressure  on  each  segment  is  evaluated. 

5.3  Node  Relocation  (R- Adaptation) 

In  order  to  reduce  element  distortion  during  complex  large  deformation  of  the  workpiece  in 
metal  forming  problems,  the  r-method  of  node  relocation  is  coupled  with  the  flexible  ALE 
description.  In  this  method,  the  motion  of  ALE  nodal  points  is  guided  by  parameters  de¬ 
termined  from  an  element  distortion  criterion,  based  on  an  inverse  mapping  algorithm  as 
suggested  in  [42]. 

For  isoparametrix  QUAD4  elements,  a  measure  of  distortion  can  be  the  included  angle 
9  between  the  adjacent  edges.  However  in  some  cases,  the  isoparametric  transformation 
resulting  from  overall  element  distortion  is  not  proper,  even  though  6  is  less  than  180°. 
This  can  be  averted  by  choosing  to  represent  distortion  through  an  inverse  mapping  as 
®^gg®sted  in  [42].  The  assumption  in  this  technique  is  that  the  number  of  iterations  required 
for  convergence  of  the  inverse  mapping  process  from  (x,y)  to  coordinates  in  equation 
(5.12)  characterizes  element  distortion. 

1  _  r/i  I  ix  —  X)tXi(l  +  \  lO'l 

V  \  [  E.  Vi  Ei  Vi  g.  \  \  4j/  -  E.  2/.(l  +  i^g.ig)  j  ^  ^  ^ 

where  Xi,yi  and  ^i,gi  are  node  coordinates  in  the  physical  and  master  domains  respectively 
for  QUAD4  elements  with  bilinear  shape  functions.  Equation  (5.12)  is  solved  iteratively  for 
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Figure  5.2:  Backward  extrusion  at  15%  deformation  (a)  without  node  relocation  (b)  with  node 
relocation. 

(^,  T]).  Node  relocation  is  enforced  with  an  ALE  node  if  the  difference  between  two  consecu¬ 
tive  absolute  values  of  ^  or  t;  exceeds  1%  after  20  iterations. 

A  weighted  averaging  formula  is  then  used  to  determine  the  new  node  location  as 

Tlk=l  Vk 

X  =  —  ^  .  ,  y  =  -=:^; — — 

^k=l  Ji^  2^k=l  Jf. 

where  n  is  the  total  number  of  surrounding  elements  sharing  the  node  to  be  relocated  and 
Xk,yk  are  the  centroidal  coordinates  of  the  neighboring  element  k.  Jk  is  the  Jacobian  in 
element  k  at  the  node  to  be  relocated,  whose  inverse  is  used  as  a  weighting  function.  This  is 
seen  to  be  more  effective  than  relocation  based  on  mean  of  centroidal  coordinates,  suggested 
in  [42],  because  it  adds  the  influence  of  the  inverse  Jacobian  to  the  relocation. 

The  effect  of  this  r-adaptation  with  ALE  nodes  are  seen  in  a  backward  extrusion  simula¬ 
tion  illustrated  in  figure  5.2.  In  this  problem  the  ALE  node  at  the  punch  edge  is  made  to 
follow  this  edge  with  progress  of  deformation.  Existence  of  Lagrangian  nodes  in  the  region 
around  the  corner  ALE  node,  gives  rise  to  element  distortion  which  are  easily  overcome  by 
this  node  relocation  scheme. 

5.4  Adaptive  Mesh  Superposition  With  S-Method 

The  multi-level  mesh  superimposition  method  termed  as  the  s-method,  proposed  by  Fish 
[21,  22],  has  the  potential  of  increasing  local  resolution  in  the  solution  by  overlaying  a  por¬ 
tion  of  a  finite  element  mesh  with  a  locally  enriched  mesh.  Such  regions  may  be  only  a  very 
small  part  of  an  element,  and  subsequent  levels  of  enrichment  may  be  added  based  on  errors 
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indicators  at  existing  levels.  This  method  is  implemented  to  predict  regions  of  localized 
deformation  in  metal  forming  simulations  by  combining  the  global-local  methodologies  in 
critical  regions  with  the  feedback  approach  of  adaptive  techniques.  Level  0  corresponds  to 
the  global  mesh  and  levels  1,2  etc.  correspond  to  the  higher  levels  of  adapted  domain.  The 
extent  of  adapted  domains  at  each  level  (1,2  ..)  decrease  with  the  increase  in  level  number. 
.4dditionally  higher  level  (1,2,..)  meshes  at  each  level  may  grow  or  diminish  with  deforma¬ 
tion.  The  level  0  mesh  in  this  study  can  be  ALE,  while  the  higher  level  adapted  meshes  are 
necessarily  Lagrangian  to  depict  material  localization.  Consequently  the  higher  level  meshes 
may  convect  through  the  level  0  mesh. 


S-adaptation  in  this  study  is  initiated  when  unacceptable  local  errors  or  steep  solution 
gradients  are  identified  in  a  existing  level  mesh.  Only  1  level  of  s-adapatation  is  explained 
in  the  ensuing  discussion.  Figure  5.3  shows  global  (level  0)  and  local  (level  1)  meshes, 
where  Cl  is  the  domain  of  the  global  mesh  with  a  boundary  F.  (Cf2  )  corresponds  to 
the  locally  enriched  region  discretized  into  a  local  mesh,  and  F*^^  is  the  interface  between 
the  global  and  local  domains  such  that  F*^^  fl  F  =  0.  The  solution  in  the  entire  domain  is 
obtained  by  superimposing  incremental  displacements  (  Au  )  from  the  global  mesh  onto  a 
local  displacement  field  (  Au^  )  defined  on  fl^  within  an  increment.  The  resulting  kinematic 
relations  in  the  incremental  analysis  are 


Au 


Au^  -|-  Au^  on  C  Cl 

Au®  on  cn,ci^  nci^  =  0 


where 

Au^  =  0  on  0,°  (5.14) 

Au^  =  0  on  F^^  (5.15) 

and 

Au®  +  Au^  =  Au^  on  F“  (5.16) 

where  F“  is  the  prescribed  displacement  boundary.  In  this  analysis,  the  displacement 
boundary  conditions  is  considered  to  be  accommodated  by  the  global  mesh,  leaving  Au^  =  0 
on  F“.  The  global  and  local  incremental  displacement  fields  are  discretized  using  standard 
C°  continuous  shape  functions 

Auf  =  N^Au^  (5.17) 

Aui  =  N^Aui^  (5.18) 

where  Auf^  and  Auf”^  are  nodal  displacement  increments.  The  compatibility  condition 
(  5.15),  is  enforced  by  the  relation 

=  0  on  (5.19) 
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5.4.1  Finite  element  formulations  for  global  local  analysis 


For  a  quasi-static  process  using  a  backward  Euler  time  integration  scheme,  the  principle  of 
virtual  work  is  written  in  the  final  grid  configuration  in  the  n-th  increment  as 


(7,j 


dSui 

dxj 


dCl 


tiSuidV 


(5.20) 


The  equation  of  virtual  work  is  solved  for  the  displacements  and  consequently  the  stresses, 
using  an  iteration  algorithm.  For  this  purpose,  it  is  convenient  to  represent  equation  (  5.20) 
as  a  scalar  valued  function  G  in  terms  of  global-local  variables  of  the  form 


Consistent  linearization  by  obtaining  the  directional  derivative  of  G  in  equation  (  5.21)  along 
incremental  displacement  vector  Au,  is  required  to  set  up  the  tangent  stiffness  matrix.  The 
quasi-Newton  family  of  update  algorithms,  such  as  BFGS  or  Broyden’s  algorithms  requiring 
a  first  approximation  to  the  stiffness,  are  employed  for  solving  the  equation  (  5.21).  These 
methods  iteratively  update  the  stiffness  matrix  to  provide  a  secant  approximation  to  the 
stiffness  matrix.  The  equations  once  assembled  has  the  following  structure 


where 


[Kf  [Kf  1  f  Au°  If  Af°  \ 
\Kf  [Kf  J  \  A-ji  /  -  \  AF^-  / 


dxi  deik  dxi 


dNf  d(7ij 
dxi  deik  dxi 


dxi  deik  dxi 


(5.22) 
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\ 


Figure  5.3:  An  example  of  the  original 


mesh  and  the  superimposed  patch  of  higher-order  element. 


dN^ 

dN^ 

+  Au^)— T — ~dQ, 

axj 

is  realized  in  equation  (  5.22). 


Level  0  mesh  is  created  using  bilinear  QUAD4  shape  functions,  while  higher  level  meshes 
use  hierarchical  shape  functions.  In  this  paper  two  levels  of  s-adaptation  is  considered  using 
upto  fourth  order  hierarchical  shape  functions.  Higher  order  hierarchical  shape  functions 
are  developed  as  a  tensor  product  of  integrals  of  Legendre  polynomials.  The  j-th  order 
hierarchical  shape  functions  in  one  dimension  are  generated  in  the  following  fashion: 

Pj-i{x)dx  (5.23) 

where  Pj  is  the  Legendre  polynomial  of  order  (j).  Detailed  discussion  on  the  development 
of  hierarchical  shape  functions  are  presented  in  [29,  22]. 

The  shape  functions  used  in  two  dimensions  can  be  separated  into  a  set  of  nodal,  side 
and  internal  shape  functions.  These  are  dependent  on  whether  the  set  corresponds  to  La- 
grangian  or  Serendipity  families.  The  space  spanned  by  the  hierarchical  shape  functions  for 


{aF^}  =  [  b,N^dn  +  /  UN^dV  -  [ 

'^^n+1  dVri^l 


\AF^}=  f  b,Nf'dn+  [  tMdT-  f 

The  hierarchical  structure  in  the  stiffness  matrix 


qL 


Hierarchical  Shape  Functions  for  s-adapted  domain 
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O  Node  with  T]  value  below  the  cutoff 

•  Node  with  T]  value  above  the  cutoff 

Figure  5.4:  Various  possibilities  in  constructing  s-adapted  subregions  in  an  element. 
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any  order  of  interpolation  is  the  same  as  that  spanned  by  Lagrangian  shape  functions  for  the 
same  order.  Two  additional  issues  are  systematically  addressed  in  the  construction  of  hier¬ 
archical  shape  functions.  They  are:  (i)  the  set  of  shape  functions  be  geometrically  isotropic, 
i.e.,  stays  the  same  irrespective  of  orientation  of  the  local  coordinate  system,  and  (ii)  com¬ 
patibility  of  displacement  modes  arising  from  odd  order  hierarchical  interpolation  functions 
belonging  to  adjacent  elements  along  a  common  edge,  are  systematically  addressed.  This 
has  been  described  in  details  in  [30].  In  this  work,  level  1  mesh  is  modeled  with  quadratic 
hierarchical  shape  functions  and  level  2  mesh  is  modeled  with  cubic  or  quartic  hierarchical 
shape  functions. 


5.4.2  Construction  of  superimposed  domains 

Generation  of  superimposed  element  patches  require  a  criteria  for  dissecting  global  elements, 
and  elements  at  subsequent  levels.  This  is  done  by  a  method  based  on  the  gradients  of 
equivalent  plastic  strains  defined  as  V  gives  a  strong  indication 

of  regions  of  localized  straining  in  large  deformation  problems.  The  following  steps  are 
sequentially  implemented  in  constructing  the  s-adapted  domain. 

•  Using  the  super  convergent  patch  recovery  technique  suggested  in  [9,  10],  gradients  of 
the  equivalent  plastic  strain  V  are  evaluated  at  element  nodal  points.  A  parameter 
T]  defined  in  equation  (5.24)  is  evaluated  at  each  node. 


V 


(V  e'^P) 


rntnimum 


(V  e^p) 


maximum 


(V 


minimum 


(5.24) 


where  (V  minimum  and  (V  e^'^)maximum  are  the  minimum  and  maximum  values  of 
the  gradient  at  all  nodal  points  in  Vt. 


•  A  cutoff  value,  rjcutoff  is  specified,  based  on  the  desired  levels.  In  general  a  value 
Vcutoff  =  0.5  is  prescribed  in  this  work.  Any  node  having  a  value  greater  than  rjcutofj 
is  overlaid  with  superimposed  elements. 

•  Every  quadrilateral  element  in  the  global  mesh  is  subdivided  into  two  triangles.  This 
eliminates  the  possibility  of  generating  polygons  with  more  than  four  sides. 

•  If  both  the  nodes  on  a  particular  element  edge  are  above  rjcutoff,  both  are  selected  as 
nodes  for  the  superimposed  element,  e.g.  edge  AD  in  figure  5.4(a).  If  neither  of  the 
nodes  on  the  element  edge  is  above  T/cuto//,  no  points  are  selected  from  that  element 
edge  (edge  BC  in  figure  5.4(a)).  If  one  of  the  nodes  is  above  rjcutoff  and  the  other  is 
below  (edge  AB  in  figure  5.4(a)),  a  linear  interpolation  is  used  to  determine  the  location 
of  the  new  node  in  the  s-domain.  The  linear  interpolation  is  used  only  when  adapting 
from  level  0  to  level  1. 


•  When  s-adapting  from  level  1  to  level  2,  the  above  procedure  is  modified  to  incorporate 
the  quadratic  variation  of  state  variables  along  element  boundaries.  In  this  case,  rj 
values  are  determined  for  midside  nodes  in  addition  to  the  corner  nodes.  If  t]  for  all 
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the  mid  side  nodes  are  greater  than  rjcutoff^  but  for  corner  nodes  are  less  than  rjcutoff, 
then  a  level  1  element  is  subdivided  into  4  quadrilaterals  as  shown  in  figure  5.4(e).  If 
T]  for  only  two  opposite  mid  side  nodes  is  greater  than  rjcutoff,  but  for  other  nodes  are 
less  than  rjcutoff,  then  the  level  1  element  is  subdivided  into  2  quadrilaterals  as  shown 
in  figure  5.4(f). 

•  If  the  mesh  generation  procedure  creates  two  adjacent  triangular  elements  sharing  a 
common  edge,  as  shown  in  figure  5.4(c),  a  single  condensed  quadrilateral  is  formed 
ensuring  that  the  included  angles  are  all  less  than  180°. 

In  metal  forming  simulations,  zones  of  strain  concentration  may  move  considerably  and  new 
zones  of  high  strains  appear  with  deformation.  Consequently  it  is  necessary  to  change  the 
extent  of  superimposed  domain  at  a  given  level,  together  with  the  addition  of  new  levels 
throughout  the  process  of  deformation. 

5.4.3  Aspects  of  Numerical  Integration 

Numerical  integration  of  the  local  terms  [K^]  and  coupling  terms  in  the  stiffness  matrix 
(equation  (  5.22))  is  performed  using  the  conventional  Gauss  quadrature  on  the  superimposed 
domain.  For  example,  corresponding  to  a  quadratic  hierarchical  element,  terms  in  [K^] 
are  approximately  of  order  4  and  consequently  3x3  Gauss  integration  rule  is  used  for 
[K^]  and  [K^]-  The  global  terms  [K^]  are  integrated  using  [B]  based  selective  reduced 
integration.  During  the  evaluation  of  stresses  in  [A^]  and  [K^]  at  integration  points  like  C 
of  the  hierarchical  element,  and  [K^]  at  integration  points  like  A  of  the  global  element  (figure 
5.5),  the  incremental  strain  contribution  from  the  superimposed  incremental  displacement 
field  and  the  global  displacement  field  are  added. 


5.5  Numerical  Examples 

Numerical  examples  using  the  techniques  developed  in  this  paper  are  divided  into  two  sets. 
The  first  two  problems  are  aimed  at  testing  and  validating  the  effectiveness  of  the  numerical 
model  with  s-adaptation.  Two  metal  forming  examples  are  then  conducted  using  the  r-s 
adapted  ALE  model.  The  material  used  for  all  examples  is  ASTM-A514  steel  with  the 
following  properties: 

Youngs  modulus  E  =  210.0  x  10®  MPa,  Poisson  ratio  v  =  0.30  and  initial  yield  stress  Co  =  6.93 
X  10®  MPa.  The  initial  void  volume  fraction  is  fo  =  0,  the  ultimate  void  volume  fraction 
is  fu  =  2/3  and  the  critical  void  volume  fraction  is  fc  —  0.15.  Plastic  strain  controlled 
void  nucleation  criteria  is  used  for  which  =  0.04,  s  =  0.1  and  =  0.3,  as  suggested 
by  Tvergaard  [32].  The  parameters  in  the  Gurson-Tvergaard  yield  function  in  equation 
(5.4)  are  taken  to  be  qi  =  1.5,  q2  =  1.0  and  q^  —  2.25.  The  uniaxial  true  stress  (a)  — 
logarithmic  strain  (e)  relation  is  expressed  using  a  power  law  as: 

e  =  ((To/E)((t/(Jo)“  (5.25) 


120 


Level  0  mesh  of  two  elements 


C 

B 


Gauss  integration  points  for  level  1  element 


Figure  5.5:  Numerical  integration  scheme, 
where  the  exponent  n  is  assigned  a  value  10. 


5.5.1  Plane  strain  tension  test  with  surface  imperfection 

In  this  example  a  plane  strain  tension  test  of  a  plate  with  surface  imperfection  is  analyzed. 
This  example  has  been  solved  using  an  extremely  refined  mesh  by  Tvergaard  [32]  and  He- 
instein  [43].  Figure  5.6  shows  the  entire  plate  and  the  region  considered  for  finite  element 
modeling.  The  plate  is  pulled  in  the  lateral  tension  by  applying  displacement  increments  to 
a  maximum  of  40%  engineering  strain.  An  imperfection  is  assumed  as  a  geometric  surface 
wave  in  the  form  of  a  cosine  function,  given  in  figure  5.6,  with  an  amplitude  ^/1  =  0.005,  h 
=  7  mm,  and  1  =  1  mm. 

The  simulation  is  executed  using  a  coarse  mesh  and  a  refined  mesh  of  495  and  1936 
QUAD4  elements  respectively.  Figure  5.7  shows  the  deformation  patterns  in  the  refined 
mesh  with  increasing  values  of  average  logarithmic  strain  e^.  It  is  observed  that  deviations 
from  uniform  straining  are  increasingly  prominent  with  increased  straining.  The  mesh  shows 
pronounced  localized  deformation  at  strain  levels  of  0.336.  Figure  5.8  shows  the  contour 
plots  of  equivalent  plastic  strain  and  void  volume  fraction  at  ta  =  0.336.  The  localized  plastic 
strain  and  void  volume  fraction  in  the  bands  are  clearly  observed  in  this  figure.  The  corre¬ 
sponding  deformation  patterns  and  contour  plots  for  the  coarse  mesh  are  depicted  in  figures 
5.9  and  5.10  respectively.  The  mesh  does  not  indicate  any  sign  of  localized  deformation, 
and  the  contour  plots  are  very  diffused  with  no  localization  of  values.  Comparison  of  figures 
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5.8  and  5.10  clearly  indicate  that  the  peak  values  of  and  /  are  grossly  under-represented 
with  the  coarse  mesh.  Mesh  sensitivity  of  the  localization  phenomenon  with  the  constitutive 
models  is  apparent  in  this  problem. 

To  overcome  limitations  of  the  coarse  mesh  in  representing  localization,  a  2  -  level  s- 
adaptation  is  activated.  The  global  mesh  (level  0)  of  495  elements  is  superimposed  with  level 
1  (grey)  and  level  2  (black)  adapted  regions,  as  illustrated  in  figures  5.11.  The  simulation 
adaptively  superimposes  the  first  level  (level  1)  of  hierarchical  quadratic  elements  at  Ca  = 
0.182,  at  which  the  parameter  77  exceeds  0.50  (equation  (5.24)).  It  is  important  that  the 
superposition  process  is  started  before  the  localization  process  begins,  failing  which  the  values 
of  localized  strains  will  be  significantly  lower.  Level  2  domain  of  quartic  hierarchical  elements 
is  superimposed  at  =  0.2231,  at  which  Tj  exceeds  0.70.  It  can  be  seen  that  level  2  region  is 
significantly  narrower  than  the  level  1  region  indicating  convergence  to  the  localized  bands. 
The  corresponding  contour  plots  at  Ca  =  0.336  are  shown  in  figure  5.12.  A  comparison  of 
figures  5.8  and  5.12  shows  that  the  results  of  coarse  mesh  with  s-adaptation  are  in  very  good 
agreement  with  those  for  the  highly  refined  mesh.  Though  the  localized  band  with  refined 
mesh  is  narrower  than  that  with  level  2  s-adaptation,  it  is  expected  that  subsequent  levels 
can  reduce  the  discrepancy.  In  figure  5.13,  a  plot  of  the  maximum  principal  logarithmic 
strain  at  a  point  A  as  a  function  of  increasing  average  strain  ta  is  shown.  Results  obtained 
with  the  refined  mesh,  coarse  mesh,  coarse  mesh  with  s-adaptation  are  compared  in  this 
plot.  The  same  results  are  also  generated  using  ABAQUS  with  the  same  refined  mesh.  It 
is  seen  that  s-adaptation  can  drastically  reduce  the  error  incurred  with  a  coarse  mesh  alone, 
especially  at  higher  values  of  ta,  thereby  producing  acceptable  results  in  comparison  with 
refined  models. 

5.5.2  Tension  test  of  a  notched  specimen 

In  this  example  a  notched  tensile  bar  under  plane  strain  uniaxial  stress  conditions  is  analyzed 
as  shown  in  figure  5.14(a).  The  specimen  is  pulled  to  a  maximum  of  6%  engineering  strain  in 
the  vertical  direction.  Only  the  upper  right  quarter  of  the  problem  is  modeled  from  symmetry 
considerations.  The  problem  is  simulated  using  a  coarse  mesh  of  450  QUAD4  elements,  a 
refined  mesh  of  1800  Quad4  elements  as  shown  in  figure  5.14(b),  and  a  coarse  mesh  with  s- 
adaptation  shown  in  figure  5.16.  Figure  5.14(b)  shows  the  refined  mesh  configuration  at  6% 
deformation,  for  which  the  contour  plots  of  equivalent  plastic  strain  and  void  volume  fraction 
are  depicted  in  figure  5.15.  Figure  5.16  shows  the  deformed  mesh  configurations  with  level 
l(grey)  and  level  2(black)  s- adaptations  at  3.2%  and  6%  strains  respectively,  corresponding 
to  77  >  0.5.  Contour  plots  are  depicted  in  figure  5.17.  From  figures  5.17  and  5.15,  it  is  seen 
that  the  peak  values  of  equivalent  plastic  strain  and  void  volume  fraction  with  s-adapted 
coarse  mesh  are  within  4%  and  10%  of  the  refined  mesh  predictions.  This  is  as  opposed  to 
24%  and  46%  difference  with  the  coarse  mesh  alone.  The  evolution  of  equivalent  stress  at 
point  A  (figure  5.14(a)  )  is  plotted  as  a  function  of  %  elongation  in  figure  5.18(a).  This  plot 
shows  that  though  the  stress  concentration  initially  increases  with  strain,  it  drops  of  quickly 
to  very  low  values  with  development  of  localized  straining  near  the  notch.  The  corresponding 
plot  (figure  5.18(b))  of  void  volume  fraction  shows  a  significant  increase  in  the  local  void 
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Figure  5.8:  Contour  plots  for  (a)  equivalent  plastic  strain  and  (b)  void  volume  fraction  at  average 
strain  of  0.336 
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Figure  5.9:  Coarse  mesh  configurations  at  average  strain  levels  of  (a)  0.216  (b)  0.336 
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Figure  5.10:  Contour  plots  for  (a)  equivalent  plastic  strain  and  (b)  void  volume  fraction  at  average 
strain  of  0.336 


Figure  5.11:  Coarse  mesh  configurations  with  s-adaptation.  (a)  initial  global  mesh  (b)  with  level  1 
of  quadratic  elements,  at  average  strain  of  0.221  (c)  with  level  2  of  quartic  elements,  and  level  1  of 
quadratic  elements,  at  an  average  strain  of  0.336 
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Figure  5.12;  Contour  plots  for  (a)  equivalent  plastic  strain  and  (b)  void  volume  fraction  at  average 
strain  of  0.336 

volume  at  A. 


5.5.3  Metal  forming  examples 

Plane  strain  analysis  is  carried  out  for  two  metal  forming  problems.  Each  of  the  problems 
are  solved  using  the  ALE  mesh  with  r-adaptation  at  three  different  resolutions  viz.  (a) 
coarse  mesh,  (b)  coarse  mesh  with  s-adaptation  and  (c)  refined  mesh.  Only  normal  contact 
is  assumed  in  these  simulations. 


Backward  Extrusion 

In  this  problem,  a  rigid  square  punch  is  made  to  back  extrude  a  rectangular  billet  inside  a 
rigid  die  as  shown  in  figure  5.19.  The  advantage  of  the  ALE  description  for  surface  nodes 
that  contact  at  the  edges  of  the  punch,  is  realized  in  this  example,  since  the  Lagrangian 
description  would  lead  to  severe  compromises  in  accuracy  [3,  5].  The  punch  is  moving  ver¬ 
tically  down  at  a  speed  of  50  cm/sec  for  0.026  seconds.  Faces  DA,  AB  and  CD  are  contact 
surfaces  on  the  rigid  punch,  EF  and  GH  are  contact  faces  on  the  die.  Only  the  right  half 
of  the  problem  is  modeled  because  of  the  symmetry.  The  effect  of  r-adaptation  in  this 
problem  is  discussed  in  section  3.  If  all  the  nodes  under  the  punch  are  ALE  nodes,  which 
have  no  relative  motion  with  respective  points  on  the  punch,  serious  mesh  distortion  can 
ensue  as  shown  in  figure  5.20.  Even  r-adaptation  scheme  cannot  avert  this  problem,  since 
elements  become  excessively  elongated.  This  is  overcome  with  by  a  special  scheme  as  follows. 
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Figure  5.13:  Maximum  principal  logarithmic  strain,  e,  at  a  material  point  A  vs  average  strain,  e, 


Figure  5.14:  (a)  Notched  tension  specimen  (b)  Refined  mesh  configuration  at  6%  elongation 


1.510E+00 


h  1.208E+00 


k  9.062E-01 


6.041  E-01 


h  3.021  E-01 


L  O.OOOE+00 


5.416E-01 


4.353E-01 


3.289E-01 


h  2.226E-01 


1.163E-01 


1.000E-02 


(a)  (b) 

Figure  5.15:  Contour  plots  for  (a)  equivalent  plastic  strain  and  (b)  void  volume  fraction  at  6.0% 
elongation  for  the  refined  mesh 
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Figure  5.16:  Coarse  mesh  configurations  with  s-adaptation  (a)  with  level  1  of  quadratic  elements 
at  3.2%  elongation  (b)  with  level  2  of  cubic  elements,  and  level  1  of  quadratic  elements,  at  6.0  % 
elongation 
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Figure  5.17:  Contour  plots  for  (a)  equivalent  plastic  strain  and  (b)  void  volume  fraction  at  6.0  % 
elongation  for  the  coarse  mesh  with  s-adaptation. 
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Figure  5.18;  Evolution  of  (a)  equivalent  stress  and  (b)  void  volume  fraction  at  A. 
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•  All  nodes  under  the  punch  are  initially  Lagrangian  (move  with  the  material),  with  the 
exception  of  an  ALE  corner  node. 

•  As  the  adjacent  Lagrangian  node  approaches  the  corner  node,  the  latter  is  released  and 
made  a  Lagrangian  node.  The  adjacent  node  is  now  moved  to  the  corner  and  made 
ALE. 

•  The  above  steps  are  repeated  throughout  the  simulation  and  r-adaptation  is  used  to 
smoothen  the  transition  between  and  Lagrangian  nodes. 

Figure  -5.21  illustrate  the  results  of  simulation  using  a  refined  mesh  with  1500  ALE 
elements.  The  corresponding  results  by  a  coarse  mesh  simulation  with  s-adaptation  are 
shown  in  figure  5.22.  Level  1  of  quadratic  superimposed  elements  (grey)  is  started  at  20% 
deformation  under  the  punch.  With  the  progress  of  deformation,  regions  of  intense  straining 
grows  and  changes,  and  consequently  the  s-adapted  level  1  region  extends.  It  is  important 
to  note  that  the  superimposed  s-domain  may  convect  through  the  level  0  ALE  mesh,  as  can 
be  seen  in  figure  5.22(b).  Level  2  mesh  superposition  begins  at  40%  deformation.  The 
superimposed  mesh  domain  in  figure  5.22(b)  gives  an  indication  of  localized  straining  under 
the  punch,  starting  at  the  corner.  Contour  plots  of  equivalent  plastic  strain  and  void  volume 
fraction  in  figure  5.22(c)  and  (d)  exhibit  good  agreement  with  refined  mesh  predictions. 
Figure  5.23(a)  show  the  contact  stress  distribution  under  the  punch  at  65%  deformation, 
while  figure  5.23(b)  presents  the  load-displacement  diagram  for  the  back  extrusion  problem. 


Ironing 

In  this  last  example,  a  rigid  prismatic  punch  is  used  to  form  a  thin  walled  hollow  cup  by 
ironing  as  shown  in  figure  5.24.  The  punch  moves  vertically  down  at  a  speed  of  50  cm/sec 
for  a  total  simulation  time  is  0.2  seconds.  Faces  AB,  BC,  CD  and  DE  are  contact  surfaces 
on  the  rigid  punch,  FG  and  HI  are  contact  faces  on  the  die  as  shown.  The  top  face  of  the 
workpiece  is  roller  supported  to  imitate  the  effect  of  friction.  Figure  5.25(a)  and  (b)  show 
the  initial  and  final  refined  mesh  (1500  QUAD4  elements)  configurations  corresponding  to 
a  punch  travel  of  7.5  cm  (75%  of  the  initial  height).  Contour  plots  of  equivalent  plastic 
strain  and  void  volume  fraction  at  the  final  stage  are  depicted  in  figures  5.25(c)  and  (d) 
respectively.  Similar  figures  for  the  coarse  mesh  simulation  (450  QUAD4  elements)  with 
s-adaptation  are  shown  in  5.26.  Level  1  (grey)  superimposition  with  quadratic  hierarchical 
elements  is  initiated  when  the  punch  travel  is  1.5  cm,  for  rj  >  0.5.  As  seen  in  figure  5.26(b), 
level  1  superimposed  domain  grows  with  deformation  in  the  form  of  a  band.  Level  2  (black) 
superimposed  cubic  hierarchical  elements  are  introduced  when  the  punch  travel  is  5  cm. 

5.6  Conclusions 

This  work  entails  two  types  of  adaptations  in  simulating  large  deformation  metal  forming 
problems.  They  are  (a)  the  r-method  of  node  relocation  for  maintaining  mesh  and  mod- 
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Figure  5.19:  Schematic  view  of  the  backward  extrusion  problem. 


Figure  5.20:  Backward  extrusion  at  40%  deformation  with  node  relocation  and  ALE  nodes  on 
contact  surface  which  are  not  released. 
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Figure  5.21:  Backward  extrusion  results  at  65%  deformation  under  the  punch  (a)  refined  mesh 
configuration  (b)  material  mesh  configuration  (c)  equivalent  plastic  strain  contours  and  (d)  void 
volume  fraction  contours 
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(a) 


(b) 

Figure  5.23:  (a)  Normal  stress  distribution  along  the  contact  region  at  65%  deformation  and  (b) 
Load-displacement  for  back  extrusion. 
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Figure  5.24:  Schematic  view  of  the  ironing  process 
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Figure  5.26:  Coarse  mesh  configuration  with  s-adaptation  (a)  at  15%  deformation  (b)  at  75% 
deformation.  Contour  plots  at  75%  deformation  for  (c)  equivalent  plastic  strain  and  (d)  void 
volume  fraction 


eling  quality,  and  (b)  a  multi-level  mesh  superposition  based  s-adaptation  for  increasing 
the  accuracy  of  local  predictions.  R-adaptation  in  conjunction  with  the  ALE  kinematic  de¬ 
scription,  is  very  effective  in  removing  element  distortions  and  better  representation  of  the 
contact  boundary.  For  metal  forming  problems,  this  helps  in  executing  the  simulation  to 
fairly  high  degrees  of  deformation  without  mesh  collapse.  The  resolution  of  a  global  finite 
element  mesh  is  increased  locally  by  s-adaptation  through  superimposition  of  hierarchical 
elements  in  regions  of  high  solution  gradients.  Strain  localization  ocurrs  over  very  narrow 
regions  in  comparison  with  the  overall  dimensions  and  the  prediction  of  localized  bands  is 
mesh  sensitive.  For  meshes  which  are  not  refined,  the  simulation  is  incapable  of  capturing 
any  strain  localization  phenomenon  at  all.  Here  the  s-adaptation  is  particularly  important 
since  localization  zones  are  apriori  unknown.  S-adaptation  in  this  study  helped  capture  lo¬ 
calized  deformations  by  superimposing  level  1  of  quadratic  hierarchical  elements  and  level 
2  of  cubic  or  quartic  hierarchical  elements.  Subsequent  levels  can  zoom  in  on  the  localized 
process.  Accuracy  of  the  solution  is  considerably  improved  with  s-adaptation  in  all  problems 
considered  in  this  study. 
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Chapter  6 


An  Arbitrary  Lagrangian-Eulerian 
Finite  Element  Model  for  Heat 
Transfer  Analysis  of  Solidification 
Processes 


Summary 

In  this  work,  a  heat  transfer  analysis  for  solidification  problems  has  been  conducted  to 
evaluate  the  temperaturefield  and  the  location  of  the  phase  change  interface.  An  arbitrary 
Lagrangian-Eulerian  kinematic  description  has  been  utilized  in  the  finite  element  formulation 
for  imparting  flexibility  to  the  motion  of  the  nodes  in  model.  By  detaching  the  nodal 
points  from  the  underlying  material,  nodes  can  be  monitored  to  follow  the  evolving  front 
while  maintaining  shapes  of  the  elements.  Special  numerical  techniques  to  smoothen  the 
deforming  front  and  to  avoid  continuous  remeshing  are  introduced.  Numerical  examples 
have  been  solved  to  establish  the  validity  of  the  present  model  and  its  strength. 


6.1  Introduction 

Accurate  modeling  of  various  solidification  processes  is  of  prime  interest  to  the  casting  in¬ 
dustry  for  enhancing  the  quality  of  Ccistings  and  reducing  the  number  of  rejects.  In  depth 
process  analysis  enables  an  engineer  to  understand  the  effects  of  temperature  history  and  gra¬ 
dients,  solid-liquid  interface  motion,  fluid  velocity  etc.  on  the  evolution  of  void  and  porosity 
distribution,  shrinkage  and  microstructure  of  solidified  body.  Design  and  control  of  various 
process  parameters  can  then  be  exercised  to  optimize  these  variables  in  order  to  obtain  the 
desired  shape,  size  and  thermophysical  properties  of  the  finished  product.  A  rigorous  thermal 
analysis  of  various  solidification  processes  entails  a  coupling  of  fluid  flow,  heat  transfer  and 
motion  of  the  phase  transition  zone.  Consideration  should  be  given  to  various  interacting 
factors  such  as; 
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a.  Transient  behavior  of  the  resulting  temperature  field. 

b.  Apriori  unknown  motion  of  the  irregular  phase  change  interface. 

c.  Effects  of  convection  in  the  liquid  pool. 

A  number  of  exact  and  approximate  analytical  studies  have  been  published  on  the  solu¬ 
tion  of  moving  boundary  problems  with  applications  in  solidification  or  melting,  and  a  cogent 
review  has  been  presented  by  Crank  [1].  Notable  among  these  methods  are  the  heat  balance 
integral  methods  by  Goodman  [2],  embedding  methods  by  Boley  [-3]  and  Greens  function 
methods  by  Hansen  and  Hougaard  [4].  Although  mathematically  elegant,  these  analyti¬ 
cal  methods  have  essentially  been  limited  to  one  dimensional  problems  with  over-simplified 
boundary  conditions  and  material  properties.  Realistic  problems  are  seldom  unidirectional 
with  simple  features,  and  hence,  cannot  be  handled  by  these  methods.  A  number  of  numeri¬ 
cal  techniques  have  been  consequently  developed  for  solving  multi-dimensional  heat  transfer 
problems  with  phase  change.  Among  them,  the  finite  difference  method  was  the  first  to  gain 
popularity  with  researchers  for  its  simplicity  in  dealing  with  partial  differential  equations.  A 
comprehensive  analysis  of  various  finite  difference  methods  that  have  evolved  in  solidification 
analysis  has  been  presented  by  Crank  [5].  Among  the  important  contributions  in  this  area 
are  the  works  of  Lazaridis  [6],  Shamsunder  and  Sparrow  [7],  Wilson  et.al.  [8]  etc.  However, 
this  method  lacks  the  generality  for  treating  practical  problems,  especially  when  complicated 
geometries  and  boundary  conditions  are  involved.  The  last  decade  has  seen  an  increasing 
trend  towards  the  development  of  the  finite  element  method  for  heat  transfer  problems  be¬ 
cause  of  its  effectiveness  with  arbitrary  shaped  domains  and  complex  boundary  conditions. 
A  classification  of  different  techniques  based  on  the  finite  element  method  is  presented  in 
Crank  [1].  A  major  challenge  in  the  analysis  of  phase  change  problems  is  the  accurate  repre¬ 
sentation  of  continuously  moving  solid-liquid  interfaces,  for  which  velocity  and  heat  transfer 
conditions  have  to  be  met.  Based  on  the  procedure  for  treating  the  release  of  latent  heat  of 
solidification  at  these  interfaces,  the  methods  have  been  traditionally  categorized  into  tw'o 
classes,  namely  the  fixed  grid  or  moving  grid  methods. 

In  the  class  of  fixed  domain  finite  element  analysis,  the  enthalpy  method  has  been  ex¬ 
tensively  used  because  of  its  relative  ease  of  implementation  in  computer  codes.  In  this 
method,  the  thermal  properties  in  the  solid  and  liquid  phases  and  the  identification  of  the 
phase  change  interface  are  all  unified  through  the  introduction  of  a  discontinuous  enthalpy 
function  H(T),  where  T  is  the  temperature  field  .  Such  an  implicit  representation  of  the 
latent  heat  release  conditions  inherently  reduces  the  burden  of  continuous  tracking  of  the 
solidification  front  and  solving  separate  equations  for  each  phase.  Among  the  large  volume 
of  published  research  on  the  enthalpy  formulation  in  finite  difference  or  finite  element  meth¬ 
ods  are  the  innovative  works  of  Crowley  [9],  Comini  et.  al.  [10],  Voller  and  Cross  [11], 
Shamsunder  and  Sparrow  [7],  Samands  et.  al.  [12]  and  Salcudean  et.  al.  [13].  Although 
the  enthalpy  finite  element  formulation  enjoys  good  convergence  rates  to  the  weak  solution 
of  the  boundary  value  problems,  they  suffer  from  some  limitations  especially  in  modeling 
isothermal  solidification  as  discussed  by  Bell  [14]  and  Voller  and  Cross  [15].  These  methods 
do  not  explicitly  account  for  the  sharp  discontinuities  and  mechanical  boundary  conditions 
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at  the  interface.  Hence,  the  accuracy  of  the  solution  deteriorates  in  the  vicinity  of  the  inter¬ 
faces.  Mesh  refinement  and  introduction  of  interfacial  control  volumes  have  been  proposed 
in  [15]  for  accuracy  enhancement,  but  have  resulted  in  rapid  decline  in  the  efficiency  of  the 
numerical  analysis. 

The  moving  grid  finite  element  methods  use  front  tracking  techniques  to  follow  the  dis¬ 
crete  phase  change  interface  continuously,  and  the  latent  heat  release  is  treated  as  a  moving 
boundary  condition.  In  these  methods,  the  energy  equations  are  written  in  terms  of  tem¬ 
perature  as  the  dependent  variable.  In  general,  it  is  accepted  that  the  moving  grid  finite 
element  methods  are  more  accurate  in  predicting  the  motion  of  the  interface  and  the  tem¬ 
perature  gradients  across  it.  Various  methods  like  deforming  grids,  variable  or  coordinate 
transformation,  introduction  of  special  algorithms  near  the  phase  change  interface  etc.  have 
been  experimented  with.  In  order  to  trace  the  free  boundary  of  the  phase  interface,  Bonnerot 
and  Jamet  [16]  have  used  a  space-time  formulation,  to  allow  for  a  change  of  nodal  positions 
with  time.  Yoo  and  Rubinsky  [17]  have  solved  a  static  heat  transfer  problem  followed  by 
internal  boundary  relocation  ,  based  on  the  solution  of  energy  balance  on  the  phase  change 
interface.  Miller  [18]  has  implemented  principles  of  error  minimization  to  move  elements  in 
Stefan  problems.  Significant  contributions  to  the  finite  element  analysis  of  moving  bound¬ 
ary  problems  have  been  made  by  Lynch,  ONeill  and  coworkers  [19,  20,  21,  22,  23].  Lynch 
and  Gray  [19]  have  developed  a  general  Galerkin  finite  element  formulation  to  automatically 
account  for  continuous  mesh  deformation  in  solving  shallow  water  wave  problems.  In  their 
papers  [20,21]  Lynch  and  ONeill  have  incorporated  implicit  time  dependence  in  the  element 
interpolation  functions  for  solving  parabolic  problems  with  phase  change.  The  finite  element 
mesh  is  allowed  to  deform  continuously  for  conforming  with  the  evolving  solid-liquid  inter¬ 
face.  Through  the  introduction  of  a  deforming/growing  numerical  coordinate  system,  the 
motion  of  nodal  points  is  explicitly  accounted  for  in  the  weak  form  of  the  governing  equations. 
No  loss  of  energy  is  therefore  encountered  while  updating  variables,  during  mesh  changes  in 
every  time  step.  This  makes  it  a  more  attractive  method  compared  to  the  static  methods 
using  remeshing  (e.g.  [17]),  in  terms  of  efficiency  and  accuracy.  Equivalence  between  vari¬ 
ous  approaches  viz.  weighted  residual,  coordinate  transformation,  space-  time  elements  [16] 
and  error  minimization  [18],  in  formulating  deforming  grid  finite  element  analysis  has  been 
established  in  a  comprehensive  paper  by  Lynch  [22].  Ensuing  shortcomings  of  the  deforming 
grid  formulation  from  excessive  mesh  distortion  were  surpassed  in  a  paper  by  Albert  and 
ONeill  [23],  where  they  monitored  the  motion  of  the  grid  by  transfinite  mapping  techniques. 

The  present  work  introduces  an  arbitrary  Lagrangian-Eulerian  (ALE)  kinematic  descrip¬ 
tion  in  the  finite  element  method  for  heat  transfer  analysis  of  problems  involving  phase 
change.  This  formulation  treats  the  finite  element  mesh  as  a  reference  frame  which  may  be 
moving  with  an  arbitrary  velocity  in  the  laboratory  system.  Thus  it  combines  the  advantages 
of  pure  Lagrangian  and  Eulerian  systems  into  one  description.  The  ALE  formulation  was 
introduced  in  the  finite  element  method  by  Donea  [24]  and  Belytschko  [25,  26]  for  modeling 
fluid-structure  interaction.  In  these  problems,  the  fluid  was  treated  with  ALE  formulation 
and  the  solid  with  pure  Lagrangian  formulation.  Liu  and  coworkers  [27]  have  extended  this 
formulation  for  three  dimensional  liquid  storage  systems.  Recently,  the  ALE  description 


146 


has  been  successfully  implemented  in  nonlinear  solid  mechanics  problems  by  Liu,  Belytschko 
and  coworkers  [28,29],  Ghosh  and  coworkers  [30,31,33]  and  Haber  [32]  for  modeling  large 
deformation.  These  numerical  models  have  been  particularly  effective  in  overcoming  some 
of  the  major  shortcomings  of  pure  Lagrangian  and  Eulerian  formulations  in  metal  forming 
simulations  [31,33],  such  as  excessive  distortion  of  elements  or  improper  representation  of 
contact  boundary  conditions. 

Though  the  arbitrary  Lagrangian-Eulerian  finite  element  method  has  some  obvious  sim¬ 
ilarities  with  the  deforming  finite  element  formulation  [22],  especially  with  respect  to  the 
development  involving  coordinate  transformation,  it  is  perceived  of  as  being  a  more  general 
and  flexible  system  that  can  be  used  for  a  wide  variety  of  applications.  Since  this  formula¬ 
tion  invokes  conservations  laws  of  mechanics  with  respect  to  an  arbitrarily  moving  reference 
system,  the  notion  of  independent  grid  motion  is  fundamentally  incorporated  into  the  finite 
element  analysis.  With  this  formulation,  different  parts  of  the  computational  domain  can  be 
represented  using  different  kinematic  descriptions.  For  example,  in  a  solidification  problem 
involving  fluid  convection,  pure  Eulerian  nodes  can  model  the  fluid  domain,  pure  Lagrangian 
nodes  can  represent  the  evolving  solid  domain  and  a  set  of  ALE  nodes  can  follow  the  mov¬ 
ing  phase  change  front  within  the  same  computational  mesh.  It  is  this  inherent  flexibility 
that  makes  the  ALE  formulation  an  attractive  candidate  for  this  class  of  problems.  In  this 
work,  a  heat  transfer  analysis  of  solidification  problems  using  this  method  has  therefore  been 
pursued. 


6.2  Governing  Equations  In  Arbitrary  Lagrangian-Eulerian  De¬ 
scription 


The  arbitrary  Lagrangian-Eulerian  description  introduces  a  reference  configuration  which 
consists  of  grid  points  in  arbitrary  spatial  motion.  Each  point  of  this  reference  configuration 
is  unambiguously  identified  by  an  invariant  set  of  three  independent  coordinates  Xi-  Concur¬ 
rently,  the  material  coordinate  system  is  represented  by  an  invariant  set  of  coordinates  JA, 
attached  to  a  material  point  and  the  spatial  coordinate  denoting  the  location  of  each  point 
in  space  is  given  by  Xi  .  The  motion  of  grid  and  material  points  can  then  be  expressed  as 
continuous  function  of  time  t  and  Xi  respectively  as  : 


Xi  =  MXj,i)  (6.1) 

=  (6.2) 

where  (f)i  and  are  mapping  functions  with  nonvanishing  Jacobians.  The  material  time 
derivative  of  any  time  dependent  physical  quantity  (3  may  then  be  related  to  the  referential 
time  derivative  by  the  relation 


dV 


=  f  I.  +  (vn  -  H 


dxk 


(6,3) 


where  Is  the  material  time  derivative  and  fracddt\^  is  the  referential  time  deriva¬ 
tive.  In  equation  (6.3)  K[=  ^\x  is  the  material  velocity  and  W[=  is  the  grid  velocity. 
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Using  the  above  relations,  material  conservation  laws  with  respect  to  an  arbitrary  reference 
system  may  be  expressed  as  : 


Continuity: 


Momentum: 


Energy: 


+  WK-vv',))  =  o 


„  ,  dcii  dVi  I  ^ 


p-§i\x  +  fK  - 


Cy jiD  ji 


where  p  refers  to  the  materiaJ  density,  dij  is  the  Cauchy  stress  tensor,  G,-  is  the  body 
force  per  unit  mass,  u  is  the  specific  internal  energy,  qi  is  the  heat  flow  vector,  Dij  is  the  rate 
of  deformation  tensor  and  /  is  the  rate  of  heat  addition  by  a  source.  Expressing  the  specific 
internal  energy  function  as  : 

U  =  C^{Q  -  Oref)  (6.7) 

where  Cp  is  the  specific  heat  at  constant  pressure,  0  is  the  absolute  temperature  and  0re j 
is  a  reference  temperature.  Incorporating  the  Fourier’s  law  of  heat  conduction,  the  energy 
balance  equation  (6.6)  may  be  restated  as  : 


+  pCpiV,-Wj)- 


d  do 


where  kij  is  the  thermal  conductivity.  The  convection  terms  containing  the  relative  velocity 
{Vi  —  Wi)  in  equations  5.4,  5.5  and  5.8,  control  the  nature  of  pointwise  description  in  the 
finite  element  model.  Thus,  when  W  is  set  to  zero,  the  description  becomes  locally  Eulerian 
while  it  is  Lagrajigian  in  the  event  that  W  =  V  at  a  point.  The  energy  balance  equation  in 
the  solid  domain  ^^(t)  is  obtained  from  equation  (6.8)  by  substituting  (Gp)«,  (A;,j)s  and  ps 
for  specific  heat,  thermal  conductivity  and  density  in  the  solid  and  that  in  the  liquid  domain 
fl;(t)  is  obtained  by  substituting  the  respective  liquid  properties  {Cp)i,{kij)i  and  pi.  The 
boundary  of  the  solid  domain  is  assumed  to  be  the  combination  of  two  basic  conditions  i.e. 
r(f)  is  ri(t)  Ur2(t)  and  may  be  delineated  as 

9{x,t)  =  6g{x,t)  y  X  on  ri(f) 


kij—ni=qn{x,t)  y  X  on  T2{t)  (6.9) 

where  corresponds  to  the  outer  unit  normal  on  the  boundary.  The  initial  domain  is 
assumed  to  be  comprised  wholly  of  the  liquid  domain  i.e. 


n,{t  =  0)  =  0  Qiit  =  0)  =  a 


(6.10) 
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where  flo  is  prescribed.  The  initial  temperature  in  the  liquid  is  assumed  to  be  specified 
everywhere  as 

0(x,O)  =  0/(x)  V  X  on  Q.i{t  =  0)  (6-11) 

The  Stefan  condition  on  the  solidification  front,  which  governs  the  a-priori  unknown 
motion  of  the  phase  change  interface  is  expressed  as 

,,  89  dd  ds 

=^i-.n  (6.12) 

where  n  denotes  the  latent  heat  of  fusion,  s(t)  refers  to  the  location  of  the  solidification 
front  and  n  is  the  unit  normal  on  the  interface  Ts-^  pointing  into  the  liquid.  Additionally, 
the  isothermal  conditions  are  assumed  on  the  solidification  front,  i.e. 

9{x,t)  =  6m  y  X  on  Ts-iit)  (6.13) 

where  9m  is  the  melting  temperature.  On  account  of  the  unknown  motion  of  the  phase 
change  interface,  equation  (6.12)  renders  the  problem  nonlinear,  even  if  the  material  proper¬ 
ties  are  temperature  independent.  Temperature  fields,  position  and  velocity  of  the  interface 
with  the  progress  of  time  are  now  the  solution  variables  of  the  problem. 


6.3  Finite  Element  Formulations  For  Heat  Transfer  Analysis 

Since  this  work  deals  with  heat  transfer  analysis  only,  the  weak  forms  of  the  mass  and 
energy  balance  equations  will  only  be  considered  here.  Such  weak  forms  are  obtained  by 
taking  the  product  of  (6.4)  and  (6.6)  with  appropriate  weighting  functions  and  integrating 
over  the  current  grid  configuration.  However,  in  the  finite  element  analysis  of  time  dependent 
problems,  emphasis  is  placed  on  the  proper  time  integration  scheme  and  hence  the  selection  of 
appropriate  configuration  for  integration.  The  commonly  used  family  of  one  step  algorithms 
using  one  parameter,  has  been  utilized  in  this  analysis.  The  stability  behavior  of  commonly 
used  algorithms,  namely  the  generalize^  rapezoidal  family  or  the  generalized  mid-point 
family  have  shown  that,  both  the  one  (q  ,  parameter  algorithms  yield  unconditional  stability 
for  a  >  I  in  the  case  of  linear  problems  [34,  36].  However,  as  pointed  out  by  Hughes  [35,  36], 
the  trapezoidal  family  gives  unconditional  stability  for  a—  1  for  nonlinear  heat  conduction 
problems  whereas  it  is  second  order  accurate  for  a  =^.  On  the  other  hand,  the  mid-point 
is  unconditionally  stable  for  oi  >  ^  and  second  order  accurate  for  a  for  nonlinear  heat 
conduction  problems.  Consequently,  this  method  of  integration  has  been  invoked  in  the 
integration  of  the  governing  equations  in  the  Lagrangian-Eulerian  formulation  .  The  weak 
form  of  the  energy  balance  equation  is  written  in  an  intermediate  grid  configuration 
as  : 


/.  +  f  pCj,{Vj  -  w,)^9dn  -  f  k 


dxj  dxi 


dft 


In 


crjiDji9dQ  +  /  f9dQ  +  f  q, 

^(*n+a)  dUdn-^-ct) 


9dT 


(6.14) 


n-f  a) 
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where  all  variables  are  at  the  intermediate  configuration  corresponding  to  time  tn+a  with  the 
position  of  a  grid  point  being  interpolated  as 


=  (1  —  a)x"  +  forO  <  a  <1 


(6.15) 


To  establish  an  unambiguous  connection  between  the  Lagrangian  and  the  referential 
system,  a  condition  that  the  material  boundary  never  leaves  the  grid  domain  may  be  imposed 
i.e. 

(V-W).n  =  0  (6.16) 

where  n  corresponds  to  the  outward  normal  on  the  boundary.  The  first  term  on  the  right 
hand  side  of  equation  (6.14)  corresponds  to  heat  generation  due  to  mechanical  coupling 
and  has  been  neglected  in  this  analysis.  The  different  variables  are  evaluated  at  tn+a  by 
interpolating  between  their  values  at  two  consequitive  equilibrium  states  at  time  and  tn+i 


on-^a  _ 


(1  -  a)^"  +  a/r 


(6.17) 


where  /9  is  a  typical  material  variable  and  a  is  the  interpolation  parameter.  To  evaluate 
the  value  of  a  material  variable  in  the  grid  configuration  at  time  tn+a,  the  increment  of  the 
variable  A^/9  at  a  fixed  grid  point  is  equated  to  that  at  a  fixed  material  point  using 

generalized  midpoint  integration  as 


AV  =  A'”/?  +  At{Wk 


y^n+a^ 


(6.18) 


Substituting  equations  (6.17)  and  (6.18)  in  equation  (6.14)  for  the  time  derivative  and 
increment  of  temperature,  the  energy  balance  equation  takes  the  form 


r  I  -  dA^d- 

/  {— pCpA^^^  +  a[pC,{Vi  -  Wi)—e  +  % 


dA^6  de 

dxj  dxj 


f 


f 


L  -  f  [pC,(Vi 


j  j  L/ o/ 2 


(6.19) 


Remark  1.  The  heat  flow  vector  on  the  second  boundary  r2(i)  may  be  assumed  to  be  a 
combination  of  flux  due  to  convection  between  the  environment  and  the  body  and  radiation 
from  the  environment.  In  the  event  that  these  conditions  or  the  temperature  are  specified 
on  the  material,  the  boundary  nodes  may  be  made  Lagrangian. 

Remark  2.  As  opposed  to  Lagrangian  formulation,  the  ALE  description  requires  contin¬ 
uous  evaluation  of  density  at  nodal  points  because  of  mass  advection  through  each  element. 
The  midstep  value  of  density  can  be  evaluated  directly  by  solving  the  weak  form  of  the 
continuity  equation  (6.4)  with  the  application  of  the  constraint  condition  (6.16). 


/  p{yi-Wi)^dn-  I 

'  ■/fifir.+a)  dXj  Jq 
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(6.20) 
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Equation  (6.20)  is  solved  for  from  which  the  density  at  an  intermediate  step  may 

be  obtained  as  : 

(6.21) 

Nevertheless,  in  this  work,  density  has  been  assumed  to  be  the  same  everywhere  in  the 
domain  and  thus  is  independent  of  temperature.  Thus  the  above  continuity  equation  has  not 
been  solved  in  the  numerical  examples  performed  here.  The  grid  domain  in  the  intermediate 
configuration  is  discretized  into  a  set  of  bilinear  isoparametric  finite  elements  with  shape 
functions  j\^(x)  are  defined  in  each  element.  Defining  the  temperature  field  and  the  weighting 
functions  in  each  element  as 


ejx.t)  =  ^,(0A-i(x) 

e{xj)  =  6,it)N,{x)  (6.22) 

the  Galerkin  matrix  form  of  the  energy  balance  ecjuation  takes  the  form 

(E  +  a  E  +  «  E  l<ab)^-'0k  =  E  ( E  +  E  ^<ab)^b  (6-23) 

E,.  gr  E. 

where  4^0,,  is  the  nodal  temperature  increment,  is  the  total  number  of  domain  ele¬ 
ments  and  El  is  the  number  of  elements  on  the  second  boundary.  In  equation  (6.23),  the 
configuration  of  the  domain  D(t„+c,)  is  not  known  apriori  and  hence  an  iterative  solution 
procedure  is  warranted. 


6.3.1  Movement  of  the  phase  change  interface 

The  motion  of  the  solid-liquid  interface  is  governed  by  equation  (6.12)  which  may  be  solved 
for  obtaining  the  normal  velocity  of  the  front  at  each  point.  However,  this  is  not  possible 
in  the  finite  element  analysis  involving  a  discretized  phase  change  boundary,  and  hence,  a 
weaker  integral  form  of  the  Stefan  condition  as  proposed  by  Lynch  [22]  and  Albert  and  O’Neill 
[2.3]  has  been_adopted  in  this  work.  Multiplying  the  Stefan  condition  (6.12)  by  a  virtual 
temperature  6  and  integrating  over  the  solid-liquid  boundary  in  the  intermediate 

configuration,  the  resulting  weak  form  becomes 


Discretizing  the  solidifying  front  into  E^-i^  two-node  linear  finite  elements  and  interpolating 
the  velocity  ^  and  virtual  temperature  B  as 


~di_  = 

0  =  B^{t)M^{x)  (6.25) 

where  yV/g,  are  linear  interpolation  functions,  the  discretized  equation  takes  the  form 


E 


ddi 


(6.26) 


In  the  ALE  formulation,  the  velocity  of  nodes  on  the  front  are  made  to  coincide  with  the 
phase  change  velocity,  given  by  Equation  (28)  specifies  only  the  normal  component  of 
velocity  of  the  moving  front  at  nodal  points,  and  hence,  the  tangential  component  is  not 
restricted.  Thus  the  tangential  components  of  nodal  velocities  may  be  specified  by  the  user 
according  to  the  requirements  of  the  problem  especially  with  respect  to  element  geometry. 
For  example,  the  tangential  component  may  be  prescribed  to  avoid  non-convex  elements  or 
element  entanglement  at  the  solid-liquid  front.  Equations  (6.23)  and  (6.26),  together  with 
the  specified  tangential  node  velocities  at  the  front  are  then  solved  simultaneously  . 


6.4  Numerical  Implementation 

In  this  section,  a  few  of  the  salient  features  of  the  solution  procedure  have  been  highlighted. 

6.4.1  Iteration  Scheme 

As  mentioned  in  section  3,  heat  transfer  analysis  of  the  solidification  process  involves  the 
iterative  solution  of  the  energy  balance  equation  (6.23)  and  Stefan  condition  (6.26)  simulta¬ 
neously  to  yield  the  resulting  temperature  field  and  the  location  of  the  phase  change  interface. 
The  quasi-Newton  family  of  iterative  solvers,  e.g.  the  Broyden  or  BFGS  update  formulae  ( 
Matthies  and  Strang  [39])  have  been 

implemented  in  this  work  for  solution  of  the  system  of  nonlinear  equations.  These  update 
methods  require  line  search  algorithms  and  are  superlinearly  convergent  to  the  solution.  The 
steps  involved  in  the  iteration  process  are  now  outlined. 

Applying  the  finite  element  approximations  in  the  discretized  domain,  a  linearized  form 
of  the  energy  balance  equation  (6.23)  may  be  written  for  the  i-th  iteration  step  as 

in+a[K](')^©  R  (6.27) 

where  A0  =  A®0'  —  A^0‘“‘  refers  to  the  correction  in  the  incremental  grid  point  tempera¬ 
ture  in  the  intermediate  configuration.  The  column  in  the  right  hand  side  of  equation 

(6.27)  corresponds  to  the  external  load  i.e.  the  right  hand  side  of  equation  (6.23),  while  the 
column  is  the  internal  load  corresponding  to  the  left  hand  side  of  (6.23). 

a)  With  an  initial  guess  for  A0  and  the  intermediate  configuration  a  starting  stiffness 

matrix  and  right  hand  vector  _'+a-A(  p(t)  jg  calculated. 

b)  An  incremental  vector,  defining  the  direction  of  the  actual  incremental  solution  vector  is 
then  obtained  as 

0'  F‘"^)  (6.28) 

c)  A  line  search  is  performed  in  the  direction  of  to  satisfy  equilibrium  by  varying  a  scalar 
parameter  s(0  <  s  <  1)  given  as 


A^0‘  =  A®0'-‘  +  s0' 


(6.29) 
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to  make  the  projection  of  the  off  balance  loads  in  the  direction  of  0'  approximately  zero. 
The  corresponding  temperature  field  in  the  mid-step  configuration  is  given  by 

0n+a  =  Q'*  (6.30) 

d)  Substituting  the  mid-step  temperature  field  0"+"  in  the  discretized  form  of  the  Stefan 
condition  (6.26)  and  adding  input  regarding  the  tangential  components  of  grid  velocities  at 
the  front,  a  global  matrix  equation  is  established.  This  equation  (6.31)  is  now  solved  for  grid 
velocities  at  the  solid-liquid  interface  as  : 

S  =  [(K^-')r‘{F*-'}'  (6.31) 

where  S  corresponds  to  the  mid-step  velocity  of  the  interface,  [(K*"')']  is  the  stiffness  ma¬ 
trix  and  corresponds  to  the  load  vector  in  the  n  +  a  th  configuration  for  the  i-th 

iteration.  It  is  to  be  noted  that  [(K*  ^)‘]  and  {F*“^}‘  should  be  updated  at  each  iteration 
step,  because  of  the  changing  values  of  the  unit  normal  n""*""  at  the  front  and  intermedi¬ 
ate  temperatures  0  .  The  location  of  the  phase  change  interface  in  now  updated  to  the 

subsequent  intermediate  configuration  by  the  relation  : 

=  S"  +  (6..32) 

The  right  hand  side  of  equation  (6.2 1),  corresponding  to  the  out  of  balance  loads,  is  now 
recalculated  based  on  the  new  v'^alues  of  temperatures,  intermediate  configuration  and  veloc¬ 
ities. 

e)  A  convergence  check  is  done  at  the  end  of  each  iteration  based  on  temperatures,  out 
of  balance  loads  and  energy.  Divergence  of  the  solution  results  when  the  out  of  balance  loads 
exceed  their  initial  value,  and  at  this  point  a  stiffness  matrix  reformation  is  performed. 

6.4.2  Evaluation  of  Normals  to  the  Interface 

Accurate  evaluation  of  normals  and  tangents  to  the  solidifying  front  is  critical  for  better  c[ual- 
ity  of  numerical  results  in  the  discrete  problem.  Non-smoothness  in  the  discretized  front  can 
lead  to  abrupt  changes  in  the  direction  of  the  normal  at  the  same  nodal  point,  corresponding 
to  two  neighboring  elements.  This  in  turn,  results  in  distinctly  different  temperature  gradi¬ 
ents  at  nodal  points  and  adversely  affects  thenodal  velocities.  A  local  smoothing  technique 
is  therefore  suggested  to  reduce  mismatches  in  the  normal  at  locations,  where  the  physical 
boundary  is  expected  to  be  smooth.  This  method  establishes  higher  order  elements  from  the 
data  available  for  two-node  line  elements,  that  are  elsewhere  used  in  the  computations. 

Consider  a  2-node  line  element  ,  corresponding  to  the  boundary  of  a  QUAD4  plane 
element.  The  coordinates  are  interpolated  by  using  a  linear  interpolation  function  in  a 
natural  coordinate  system  (  — 1  <  s  <  1)  as: 

Xh(s)  =  x(-l)i(l  -s)  +  x(l)^(l  +  s)  (6..33) 
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Using  a  Taylor  series  expansion  in  the  neighborhood  of  x(s)  and  retaining  upto  quadratic 
terms,  the  coordinates  of  the  interface  points  near  x(s)  may  be  approximated  as  : 


x(s°)  =  x(s)  +  (s°  -  5)^(5°  -  s)  +  -  5)^1^ 


(6.34) 


Substituting  equation  (6.33)  in  (6.34)  yields  an  interpolation,  that  accounts  for  quadratic 
terms  and  is  smoother  than  (6.33). 


1 

x(s)  =  Xh(s)  +  -{s^  -  1)  — 


(6..35) 


The  second  derivative  in  equation  (37)  is  evaluated  by  using  a  method  of  weighted  averages 
from  values  in  neighboring  elements.  Due  to  its  quadratic  nature,  equation  (37)  can  be 
considered  to  represent  the  3-node  quadrilateral  line  elements  and  the  coordinates  of  a  third 
point  at  s  =  0  can  be  computed  from  it.  The  slope  to  the  interface  and  hence  the  normal 
vectors  are  then  computed  at  any  point  from  interpolation  based  on  the  three  nodal  points 
as: 


ds 


ds 


(6.36) 


6.4.3  Introduction  of  a  Pseudo  Domain  for  Discretization 

A  solidification  problem  involves  an  internal  solid-liquid  boundary  which  is  continuously 
evolving  with  time.  In  the  ALE  formulation,  a  set  of  nodes  is  placed  on  the  moving  front  at 
the  start  of  simulation  and  is  made  to  adhere  to  it  throughout  the  process.  This  calls  for  a 
continuous  rediscretization  of  the  computational  domain,  resulting  in  changing  the  degrees 
of  freedom  and  element  connectivity  of  the  finite  element  model.  For  example,  if  the  domain 
initially  consists  of  the  liquid  only,  elements  should  be  added  as  the  solidifying  front  forms 
and  moves  through  the  liquid.  To  avoid  the  burden  of  having  to  update  the  degrees  of  free¬ 
dom,  element  connectivity  and  hence  the  order  of  the  stiffness  matrix  and  load  vectors  ,  a 
pseudo-domain  is  introduced  adjacent  to  the  solidifying  front  as  shown  in  figure  1.  This  do¬ 
main  is  discretized  along  with  the  physical  domain,  to  yield  a  combined  number  of  elements 
which  remain  unchanged  throughout  the  simulation.  A  few  special  assumptions  are  imposed 
on  this  part  of  the  finite  element  model  as  follows. 

i)  The  specific  heat  capacity  Cp  is  zero,  i.e.  this  domain  is  always  in  a  steady  state  condition. 

ii)  If  the  boundary  Tactual  is  of  the  second  type,  i.e.  the  heat  flux  normal  to  the  boundary  is 
specified,  then  a  heat  flux  should  be  specified  on  the  pseudo  boundary  T pseudo  ^  This  is  done 
according  to  the  assumption  that  the  normal  heat  flux  through  the  pseudo  boundary  T pseudo 
at  any  time  t,  is  equal  to  the  heat  flux  through  the  adjacent  physical  boundary,  or 

I  qifiidT  =  I  qiTiidT  (6.37) 

The  other  sides  of  the  pseudo-domain  are  insulated. 
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Element  in  Pseudo  Domain 


Figure  6.1:  Introduction  of  the  pseudo  domain 

iii)  If  temperature  conditions  are  prescribed  on  the  actual  boundary,  from  which  solidification 

initiates,  additional  efforts  are  needed.  First  of  all,  the  prescribed  boundary  temperatures 
are  imposed  on  the  common  boundary.  Secondly,  heat  flux  at  the  nodes  of  the  common 
boundary  are  calculated  from  adjacent  element  temperatures  in  the  physical  domain.  This 
heat  flow  vector  is  then  prescribed  on  the  pseudo-boundarv  according  to  equation 

(6,37). 

iv)  The  boundary  of  the  pseudo-domain  coincides  with  the  that  of  the  physical  domain  ap¬ 
proximately  along  the  part  from  which  solidification  takes  place. 

v)  Nodes  from  the  pseudo-domain  continuously  move  into  the  physical  domain  by  virtue  of 
the  ALE  formulation  to  maintain  well  behaved  elements  in  the  physical  domain.  Since  the 
pseudo-domain  is  merely  a  node  bank,  element  shapes  in  this  region  is  not  important. 

6,5  Numerical  Examples 

A  two-dimensional  program  has  been  developed  with  QL^AD4  domain  elements  and  2-  node 
line  elements  on  the  solid-liquid  interface.  The  value  of  a  for  intermediate  configuration  has 
been  chosen  to  be  1/2.  The  results  of  the  program  have  been  compared  with  all  the  heat 
transfer  results  presented  in  Zabaras  [38].  Finally  a  simplified  form  of  a  continuous  casting 
process  has  been  simulated.  In  the  numerical  examples,  it  has  been  assumed  that  the  liquid 
domain  is  maintained  at  a  uniform  temperature  corresponding  to  the  melting  temperature 
and  there  are  no  convection  effects  in  the  liquid.  Consequently,  only  the  solid  domain  has 
been  modeled. 

In  accordance  with  the  data  presented  in  [38],  the  following  thermal  properties  have  been 
assumed  for  the  material  (Aluminum,  in  this  case  )  for  all  the  problems.  They  are 
I\s  —  0.0548  kcal/ms.  (Conductivity  of  solid) 

Ki  =  0.0548  kcal/ms  .""C  (Conductivity  of  liquid) 
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(C'p)^  =  0.2526  kcal/kg  .  °  C  (Specific  heat  of  solid) 

{Cp)i  =  0.2526  kcal/kg  .  °  C  (Specific  heat  of  liquid) 

L  =  94.44  kcal/kg  (Latent  heat  of  fusion) 
p  =  2650  kg/m  (Density  of  solid/liquid) 

0[=  660  °C  (Initial  liquid  temperature) 

0m=  660  '^C  (Melting  temperature) 

Furthermore,  as  suggested  in  [.38],  a  temperature  condition  has  been  specified  at  the  cooling 
boundary  as: 

doit)  =  ea{t)  +  {emit)0a{t))e-’^^  (6.38) 

where  9a{t)  corresponds  to  the  ambient  temperature  and  is  a  cooling  parameter. 

6.5.1  An  one-dimensional  problem 

A  one  dimensional  solidification  problem,  for  which  the  initial  domain  and  discretization  is 
shown  in  figure  2(a),  is  solved  with  the  boundary  conditions  specified.  In  equation  (6.38), 
R  is  assumed  to  have  a  value  of  0.023  and  two  different  cases  have  been  solved  with  the 
ambient  temperature  0a{t)  taking  a  value  of  300  °C  and  51.7  °C  respectively.  The  actual 
domain  is  initially  made  up  entirely  of  liquid,  and  102  QUAD4  elements  are  packed  into  the 
pseudo-domain  as  shown.  As  solidification  begins,  the  nodes  on  the  right  hand  face  of  the 
model  move  with  the  front.  The  set  of  nodes  behind  the  front  are  moved  to  replace  the  former 
set,  so  as  to  retain  the  actual  boundary  with  temperature  boundary  conditions.  Following 
this,  ALE  nodes  are  admitted  periodically  into  the  physical  domain  and  moved  to  avoid 
elongated  elements.  Figure  2(b)  shows  the  discrete  solidified  domain  and  the  temperature 
distribution  for  6a{t)  =  300  “C.  The  motion  of  the  solid-liquid  interface  as  a  function  of 
time  has  been  plotted  in  figure  3  (a).  The  solidification  front  has  a  greater  velocity  for 
the  lower  ambient  temperature,  as  can  be  physically  predicted.  Figures  3(b)  and  4  are 
plots  of  temperatures  at  three  different  locations  for  the  two  different  ambient  temperatures. 
Excellent  agreement  is  obtained  for  all  these  results  with  those  predicted  by  Zabaras  [38]. 

6.5.2  Solidification  of  a  Cylindrical  Body 

Solidification  of  a  cylindrical  body  is  considered  in  this  two  dimensional  example.  An  ambient 
temperature  0a{t)  =  500  °C  is  prescribed  along  the  circumference  of  a  circular  section  of 
radius  18  mm  and  R  =  0.1  The  simulation  is  executed  for  8.32  seconds  at  which 

time  the  solidification  is  almost  complete.  Figure  5(a)  shows  the  initial  discretization  of 
the  pseudo-domain  with  240  elements.  Only  a  quarter  of  the  problem  has  been  considered 
because  of  symmetry  with  appropriate  boundary  conditions  indicated  on  figure  5(a).  Figures 
5(b)  and  5(c)  show  the  mesh  domain  and  the  temperature  in  the  solidifying  medium,  at  times 
4.0  secs  and  8.32  secs  respectively.  The  total  number  of  elements  in  the  pseudo  and  physical 
domains  remain  the  same  at  all  times.  The  results  of  this  problem  are  plotted  and  compared 
with  those  given  in  [38].  Figure  6  (a)  is  a  plot  of  the  front  movement  as  a  function  of  time 
while  in  figure  6(b),  the  temperature  histories  at  different  radii  are  depicted.  Except  for  a 
slight  deviation  near  the  completion  of  solidification,  the  front  radius  in  figure  6(a)  matches 
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Figure  6.2:  (a)  Initial  domain  (b)  final  domain  at  time  t=  300  s  6a{t)  =  300  '"C 
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Figure  6.3:  (a)  Position  of  the  front  for  =  .51.7  °C  and  0^(1)  =  300  °C 
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Max.' 


Figure  6.0.  (a)  Initial  discretization  of  the  pseudo-domain  (b)  mesh  configuration  and  temperature 
distribution  for  time  t=4s,  and  (c)  t=  8.32  s  in  cylindrical  problem 

with  Zabaras'  results  very  well.  Also  the  temperature  profiles  are  in  good  agreement  at 
different  radial  locations. 

6.5.3  Solidification  of  a  Square  Region 

A  40mm  x  40mm  square  liquid  domain  with  a  specified  ambient  temperature  9a{t)  —  500 
on  its  boundary  is  now  simulated  for  solidification.  As  before,  only  a  quarter  of  the  square 
region  is  analyzed  because  of  symmetry  and  appropriate  boundary  conditions  are  used.  As 
in  the  previous  example,  the  value  of  R  is  chosen  to  be  0.1  The  process  is  simulated  for 
10.4  seconds  for  solidification  to  be  nearly  complete.  Figure  7(a)  shows  the  discretization  of 
the  pseudo-domain  into  166  QUAD4  elements  with  specified  boundary  conditions.  Figures 
7(b)  and  7(c)  show  the  combined  mesh  configurations  and  temperature  distributions  in  the 
solidifying  medium  at  5.0  secs  and  10.4  secs  respectively.  It  may  be  observed,  that  towards 
the  end  of  simulation,  almost  all  the  elements  in  the  pseudo-domain  have  moved  to  the 
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(a)  Radial  location  of  the  front  versus  time  (b)  temperature  history  at  various  radial 
cylindrical  solidification 
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Figure  6./.  (a)  Initial  discretization  of  the  pseudo-domain  (b)  mesh  configuration  and  temperature 
distribution  for  time  t=5  s,  and  (c)  t—  10.4  s  in  square  problem 

physical  domain.  Figure  8(a)  shows  a  plot  of  the  x-  and  y-  coordinates  of  the  front  at 
various  times  in  the  solidification  process.  In  figure  8(b),  the  temperature  histories  at  two 
different  locations  are  found  to  agree  well  with  the  results  of  Zabaras  [38]. 

6.5.4  Simplified  Example  of  a  Continuous  Casting  Process 

In  this  final  example,  a  simplified  model  for  two-dimensional  continuous  strip  casting  has 
been  simulated.  As  illustrated  in  figure  9(a),  the  liquid  is  assumed  to  move  from  left  to  right 
with  a  horizontal  velocity  of  10  mm/sec.  The  top  surface  of  the  liquid  is  in  contact  with 
water  and  accordingly  the  heat  flux  across  this  boundary  is  assumed  to  be  given  as 

,  d6 

=  -ko{9-e^ater)  (6.39) 

where  is  a  convection  coefficient  =  5  kcal/  m^  sec  "C  and  -  30  "C  .  The  right  surface 

representing  the  free  end  of  the  solid  is  assumed  to  be  in  contact  with  a  medium  with  very 
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Figure  6.8:  (a)  Location  of  the  front  versus  time  (b)  temperature  history  at  various  locations  in 
square  solidification 
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low  convection  coefficient  and  hence  the  heat  flux  across  this  boundary  is  set  to  zero.  The 
convection  boundary  condition  (6.39)  is  directly  applied  to  the  actual  boundary  at  all  times. 
The  nodes  that  are  completely  within  the  pseudo-domain  are  apriori  set  and  maintained  at 
a  specified  temperature.  The  center-line  or  line  of  symmetry  is  at  a  distance  1  mm  from  the 
top  surface  of  the  liquid.  Various  steps  of  the  solidification  are  depicted  in  figure  10  .  The 
pseudo-domain  has  been  excluded  from  these  figures  for  better  understanding  of  the  actual 
domain.  Nodes  at  the  top  boundary  of  the  liquid  domain  are  at  rest  in  each  time  step,  i.e. 
these  are  Eulerian  nodes  .  However,  as  the  front  progresses,  nodes  from  the  pseudo-domain 
move  down  to  replace  the  previous  set  of  nodes  on  this  boundary.  The  latter  set  of  nodes 
now  move  into  the  solidified  region.  All  nodes  in  the  solidified  region  are  ALE  nodes  and 
have  the  flexibility  to  be  moved  independent  of  the  material.  The  horizontal  component  of 
velocity  i  for  nodes  on  the  right  hand  boundary  of  the  solidified  region  are  maintained 
at  the  material  velocity.  The  vertical  component  of  velocities  for  ALE  nodes  in  the 
solidified  domain  including  the  solid-liquid  boundary,  are  adjusted  so  as  to  get  a  uniform 
spacing,  shown  in  figures  17  and  18.  On  the  interface,  the  horizontal  node  velocities  are 
determined  from  the  normal  velocity  IT'„,  horizontal  velocity  Wj,  and  the  direction  cosines 
of  the  local  tangent  .  When  the  front  reaches  the  center-line,  the  vertical  component  of 
velocities  are  set  to  zero.  Also  the  temperature  boundary  conditions  are  changed  to  zero 
heat  flux  conditions.  Figure  10  shows  the  mesh  movement  and  temperature  contours  at 
various  stages  of  simulation  upto  0.56  secs,  during  which  time  the  solidified  region  has  a 
length  of  5.6  mm.  Figure  9(.b)  depicts  the  x-  and  y-  coordinates  of  the  front  as  a  function  of 
time. 

6.6  Conclusions 

In  this  work  an  arbitrary  Lagrangian-Eulerian  finite  element  method  has  been  implemented 
for  solving  the  heat  transfer  problem  in  simulation  of  solidification  processes.  The  flexibilitv 
imparted  by  introducing  an  arbitrary  reference  frame,  is  favorably  used  in  monitoring  the 
solid-fluid  interface  as  well  as  in  maintaining  well  behaved  elements  in  the  solidified  region. 
Through  the  introduction  of  a  pseudo-domain  ,  continuous  addition  of  nodes  to  the  solidifving 
body  has  been  avoided.  Comparison  of  the  results  with  some  of  the  previously  reported 
results  in  [38]  shows  excellent  agreement.  A  simplified  form  of  continuous  casting  that  was 
solved  with  this  algorithm  yields  expected  nature  of  results.  Overall  the  ALE  finite  element 
method  shows  promise  for  this  class  of  problems  and  is  expected  to  be  used  for  realistic 
simulations  in  future. 
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Figure  6.9:  (a)  Initial  discretizationwith  boundary  conditions  (b)  progression  of  the  front  for  various 
times  in  continuous  Cctsting 
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